set.seed(12345) # for reproducibility
options(knitr.kable.NA = '')

# install // load packages 
# Some packages need to be loaded. 
# We use `pacman` as a package manager, which takes care of the other packages. 
if (!require("devtools", quietly = TRUE)) install.packages("devtools")
if (!require("papaja", quietly = TRUE)) devtools::install_github("crsh/papaja")
if (!require("patchwork", quietly = TRUE)) devtools::install_github("thomasp85/patchwork")
if (!require("klippy", quietly = TRUE)) devtools::install_github("RLesur/klippy")
if (!require("pacman", quietly = TRUE)) install.packages("pacman")
if (!require("Rmisc", quietly = TRUE)) install.packages("Rmisc")
if (!require("ggbeeswarm", quietly = TRUE)) install.packages("ggbeeswarm") # Never load it directly.
pacman::p_load(tidyverse, papaja, knitr, dplyr, car, psych, afex, lme4, lmerTest, 
               emmeans, ggplot2, ggpubr, lattice, latticeExtra, parallel, effects, psycho, caret)
library("patchwork"); library("klippy")
klippy::klippy()

1 Abstract


이 분석은 Repetition Suppression / Enhancement의 신경 기제를 Voxel-wise Background Connectivity를 통해 살펴보는 것에 초점을 맞추고 있다. 구체적으로, Repetition Suppression / Enhancement를 설명하는 서로 다른 두 가지 모델. Sharpening Model과 Predictive Coding Model을 비교하고자 한다.


Repetition suppression과 enhancement * 반복되는 대상에 대한 신경 활성화의 감소(Suppression) 및 증가(Enhancement) * 일반적으로 대상 영역의 전체 voxel 중 RS voxel은 75%, RE voxel은 25%의 비율로 나타난다. * 친숙하지 않은 대상, 알아보기 힘든 대상한테는 RS보다 RE가 발생, 또는 RS의 효과 감소하는 모습을 보인다.


Sharpening Account * 반복은 대상의 표상에 덜 관련성이 있는 neural population의 반응을 감소시킨다 * 대상을 처음 표상할때는 대상에 관련이 있는, tuning된 neural population 뿐만 아니라 그렇지 않은 neural population도 함께 반응하게 된다. * repetition이 발생하면 tuning되지 않은 neural population들의 반응이 감소하면서 전체적으로 RS가 발생한다. * RS는 Sharpening에 의해 selectivity가 적은 neural population의 반응 감소를 반영하며 RE는 미처 Sharpening되지 않은 population이나 대상을 표상하는데 selectivity를 가지고 있는 neural population의 반응 증가를 반영한다. * 이에 따르면, RS를 보이는 voxel들은 RE를 보이는 voxel들 보다 대상을 표상하는데 관련성이 적은 voxel들일 수 있다. 즉, RS보다 RE가 중요하다.


Predictive Coding Account * 우리의 뇌는 감각 입력을 예측하기 위해 기억에 기반한 맥락과 정보로 외부 환경에 대한 모델을 생성한다. * 상위 영역에서는 기억에 기반하여 Prediction을 생성하고, Feedback을 통해 하위 감각 영역으로 전달된다. * 반대로 하위 영역에서는 Prediction과 실제 Sensory input을 비교하여 둘 간의 Mismatch를 Prediction error signal로 부호화하고, 상위 영역으로 Feedforward한다. 이러한 상호작용을 통해 모델을 지속적으로 갱신할 수 있다. * Repetition은 대상에 대한 Prediction을 형성하는데 영향을 준다. Repetition이 더 많이 일어날수록, 상위 영역에서 더 정교한 Prediction이 형성되며 하위 영역에서의 Error signal이 감소한다. * RS는 정교한 Prediction에 의해 Error signal이 감소되는 것을 반영한다. RE는 Prediction에 대한 Error signal의 증가를 반영한다. * 다른 측면에서 보면, RS는 Prediction을 주로 받는 population, RE는 Prediction Error Signal을 주로 생성하는 population일 수 있다. * 이에 따르면, RS와 RE는 서로 다른 functional role을 가진다. Sharpening에서 제안하는 것처럼 RS가 대상을 표상하는데 관련성이 적은 voxel들인 것은 아니다. 오히려 대상을 표상하는 중요한 representation unit일 가능성이 크다.


  • RS voxel과 RE voxel 중 어떤 것이 대상을 표상하는 데 더 중요한 정보를 많이 가지고 있는가? 어떤 모델이 더 적절할까?
  • Sharpening 모델에 따르면 RS는 불필요한 신경 신호를 반영하므로 RS 보다 RE가 대상의 표상에 더 중요한 정보를 가지고 있어야 한다.
  • 반면, PC 모델에 따르면 RS는 Prediction, Representation, RE는 Error signal을 반명하므로 RE 보다 RS가 중요한 정보를 가지고 있어야 한다.


  • 이 질문은 두 가지 방식으로 검증할 수 있다.
  1. Deep learning 기법을 사용. RS voxel과 RE voxel을 나누어 얼굴 표상 재구성에 중요도를 평가한다.
  2. Voxel-wise Background Connectivity를 사용. Repetition Suppression의 효과와 상위 영역과의 연결성 간 상호작용을 검증한다. 추가로 RS voxel/RE voxel과 상위 영역 간의 연결성을 평가한다.


  • 이 분석에서는 2번째 방법을 사용한다. 가설은 다음과 같다.
  1. Sharpening. 대상을 표상하는 데 RS보다 RE가 중요. 중요한 voxel들은 상위 영역과 연결성이 클 가능성이 있음. RS는 불필요하므로 상위 영역과 연결성이 낮은데, RE는 중요하므로 상위 영역과 연결성이 높아야함.
  2. PC. RS는 Prediction & Representation, RE는 Error Signal을 반영하므로 RS가 대상 표상에 중요. RS는 상위 영역과 연결성이 높은 반면, RE는 연결성이 낮을 것임.




2 Results 1 - Repetition Suppression Effect & bgFC


2.1 FFA/PPA - VC Background Connectivity (Phase 1)


Phase 1에서 추출한 FFA ROI과 fBGRED 참가자들의 좌표로 만든 PPA ROI가 제대로 추출되었는지 확인하기 위해 Phase 1에서 FFA/PPA-VC Background Conectivity를 비교하였다. Phase 1에서 참가자들이 본 이미지 자극들은 모두 얼굴이므로, ROI Localizing이 잘 되었다면, PPA-VC bgFC 보다 FFA-VC bgFC가 더 강하게 나타나야한다.




t2 <- read.csv("data/bg_summary_nsm.csv", header = T)
t2_l <- gather(t2, hRoi, fc, ppa:ffa, factor_key=TRUE)
t2_l$subj = factor(t2_l$subj, levels=c("primeRecon_2", "primeRecon_4", 
                                       "primeRecon_5", "primeRecon_6",
                                       "primeRecon_7", "primeRecon_8", 
                                       "primeRecon_9", "primeRecon_10",
                                       "primeRecon_11", "primeRecon_13",
                                       "primeRecon_14", "primeRecon_15",
                                       "primeRecon_16", "primeRecon_17",
                                       "primeRecon_18", "primeRecon_19",
                                       "primeRecon_20", "primeRecon_21",
                                       "primeRecon_22", "primeRecon_23",
                                       "primeRecon_24"), 
                   labels=c(2,4,5,6,7,8,9,10,11,13,14,15,16,17,18,19,20,21,22,23,24))


t2_l <- t2_l %>% filter(subj!=99)
t2_l$vRoi = factor(t2_l$vROIs, levels=c("v1","v2","v3","v4","all"),
                   labels=c("v1","v2","v3","v4","v_all")) 
t2_l$hRoi = factor(t2_l$hRoi, levels=c("ppa","ffa")) 
t2_l <- t2_l %>% select(subj,hRoi,vRoi,fc)

# number of subject
length(unique(t2$subj))
## [1] 20

# variables
glimpse(t2_l, width = 70)
## Rows: 200
## Columns: 4
## $ subj <fct> 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6…
## $ hRoi <fct> ppa, ppa, ppa, ppa, ppa, ppa, ppa, ppa, ppa, ppa, ppa, …
## $ vRoi <fct> v1, v2, v3, v4, v_all, v1, v2, v3, v4, v_all, v1, v2, v…
## $ fc   <dbl> -0.281099400, -0.231826700, -0.137576500, -0.077168000,…


각 시각 영역과 상위 영역 간 Background Connectivity 결과의 요약치는 아래와 같다.



# subject-level, long format
t2_allL <- t2_l %>% group_by(subj, vRoi, hRoi) %>%
  dplyr::summarise(fc=mean(fc)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'vRoi'. You can override using the `.groups` argument.
# t2_allL %>% kable(digits=2)

# subject-level, wide format
t2_allW <- t2_allL %>% spread(key=hRoi, value = fc)
# t2_allW %>% kable(digits=2)

# summary table: grand mean
t2_allG <- t2_allL %>% group_by(vRoi, hRoi) %>%
  summarise(fc.m = mean(fc), fc.sd = sd(fc)) %>%
  ungroup()
## `summarise()` has grouped output by 'vRoi'. You can override using the `.groups` argument.
t2_allG$fc.se <- Rmisc::summarySEwithin(data = t2_allL, measurevar = "fc", 
                                        idvar = "subj", withinvars = c("vRoi", "hRoi"))$se
t2_allG$fc.ci <- Rmisc::summarySEwithin(data = t2_allL, measurevar = "fc", 
                                        idvar = "subj", withinvars = c("vRoi", "hRoi"))$ci
t2_allG <- t2_allG %>% 
  mutate(lower.ci = fc.m-fc.ci,
         upper.ci = fc.m+fc.ci,
         lower.se = fc.m-fc.se,
         upper.se = fc.m+fc.se)
# for between-subject design. check help(summarySE)
t2_allG<-t2_allG %>% select(vRoi, hRoi, fc.m, fc.sd, fc.se, fc.ci, lower.ci, upper.ci, lower.se, upper.se)
t2_allG %>% select(vRoi, hRoi, fc.m) %>% spread(hRoi, fc.m) %>% kable(digits=2)
vRoi ppa ffa
v1 -0.04 0.35
v2 0.00 0.43
v3 0.13 0.47
v4 0.19 0.46
v_all 0.04 0.47


아래 그래프는 각 시각 영역과 상위 영역 간의 Background Connectivity (Zr)를 나타낸다. 점은 평균, error bar는 95% ci를 나타내었다.


t2.all.plot2 <- ggplot(data=t2_allL, aes(x=vRoi, y=fc, fill=hRoi)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  
  geom_segment(data=filter(t2_allW, vRoi == "v1"), inherit.aes = FALSE,
               aes(x=0.80, y=filter(t2_allW, vRoi == "v1")$ppa,
                   xend=1.20, yend=filter(t2_allW, vRoi == "v1")$ffa),
               color="gray80") +
  geom_segment(data=filter(t2_allW, vRoi == "v2"), inherit.aes = FALSE,
               aes(x=1.80, y=filter(t2_allW, vRoi == "v2")$ppa,
                   xend=2.20, yend=filter(t2_allW, vRoi == "v2")$ffa),
               color="gray80") +
  geom_segment(data=filter(t2_allW, vRoi == "v3"), inherit.aes = FALSE,
               aes(x=2.80, y=filter(t2_allW, vRoi == "v3")$ppa,
                   xend=3.20, yend=filter(t2_allW, vRoi == "v3")$ffa),
               color="gray80") +
  geom_segment(data=filter(t2_allW, vRoi == "v4"), inherit.aes = FALSE,
               aes(x=3.80, y=filter(t2_allW, vRoi == "v4")$ppa,
                   xend=4.20, yend=filter(t2_allW, vRoi == "v4")$ffa),
               color="gray80") +
  geom_segment(data=filter(t2_allW, vRoi == "v_all"), inherit.aes = FALSE,
               aes(x=4.80, y=filter(t2_allW, vRoi == "v_all")$ppa,
                   xend=5.20, yend=filter(t2_allW, vRoi == "v_all")$ffa),
               color="gray80") +
  geom_point(data=t2_allL, aes(x=vRoi, y=fc, fill=hRoi), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t2_allG, aes(x = vRoi, y=fc.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("V1","V3","V3","V4","all VC")) +
  scale_fill_manual(values = c("#feb24c", "#91bfdb"), # c("#feb24c", "#91bfdb"),
                    labels = c("PPA", "FFA")) +
  coord_cartesian(ylim = c(-0.3, 0.7), clip = "on") +
  labs(x = "", y = "Background Connectivity (Zr)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
#ggsave("bg_plot1.jpg", plot = t2.all.plot2, width=8, height=8, unit='in', dpi=600)
t2.all.plot2


각 시각 영역(V1 ~ V4)과 관심 상위 영역(PPA, FFA)을 요인으로 4 x 2 혼합 요인 분산분석을 수행하였다.


t2_allL.1 <- t2_allL %>% filter(vRoi!="v_all")
t2.aov1 <- aov_ez(id="subj", dv="fc", data = t2_allL.1, within = c("vRoi", "hRoi"))
# summary(t2.aov1)
nice(t2.aov1, es="pes") %>% kable(digits=3)
Effect df MSE F pes p.value
vRoi 1.61, 30.54 0.02 28.73 *** .602 <.001
hRoi 1, 19 0.03 156.68 *** .892 <.001
vRoi:hRoi 1.31, 24.94 0.01 13.41 *** .414 <.001


모든 효과가 유의하였다.


# post-hoc - difference between hROI, by attention condition & vROI
t2.aov1.m1 <- emmeans(t2.aov1, pairwise ~ hRoi | vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t2.aov1.m1$contrasts %>% kable(digits=3)
contrast vRoi estimate SE df t.ratio p.value
ppa - ffa v1 -0.397 0.037 19 -10.796 0
ppa - ffa v2 -0.421 0.030 19 -14.008 0
ppa - ffa v3 -0.342 0.030 19 -11.452 0
ppa - ffa v4 -0.277 0.033 19 -8.486 0
# plot(t2.aov1.m1, horizontal = T, comparison=T)

# post-hoc - difference between attention, by hROI condition & vROI
t2.aov1.m2 <- emmeans(t2.aov1, pairwise ~ vRoi | hRoi , adjust = "tukey") # adjust="bon" , "tukey", "none"
t2.aov1.m2$contrasts %>% kable(digits=3)
contrast hRoi estimate SE df t.ratio p.value
v1 - v2 ppa -0.049 0.017 19 -2.858 0.046
v1 - v3 ppa -0.169 0.024 19 -6.961 0.000
v1 - v4 ppa -0.229 0.039 19 -5.797 0.000
v2 - v3 ppa -0.120 0.015 19 -7.979 0.000
v2 - v4 ppa -0.180 0.030 19 -5.986 0.000
v3 - v4 ppa -0.060 0.021 19 -2.813 0.050
v1 - v2 ffa -0.073 0.018 19 -3.965 0.004
v1 - v3 ffa -0.115 0.021 19 -5.441 0.000
v1 - v4 ffa -0.109 0.031 19 -3.505 0.012
v2 - v3 ffa -0.042 0.011 19 -3.777 0.006
v2 - v4 ffa -0.036 0.024 19 -1.490 0.462
v3 - v4 ffa 0.006 0.018 19 0.312 0.989
# plot(t2.aov1.m2, horizontal = T, comparison=T)


사후 검정 결과, 모든 시각 영역에서 PPA 보다 FFA 와 bgFC가 더 강하게 나타났다.


2.2 관심 영역의 Repetition Suppression/Enhancement


전체 관심 영역 (FFA/PPA, V1~V4) 에서 Repetition Suppression / Enhancement 효과를 검증하였다. 제시 조건은 Initial과 Repeated로 구분되었다. 이 분석에서는 전처리 과정에서 FWHM 0mm로 Spatial Smoothing이 수행되지 않은 데이터를 사용하였다.



# higher ROIs
t3 <- read.csv("data/rse_summary_hroi_nsm.csv", header = T)
t3_l <- t3
# change class of main factors: double to factor
t3_l$subj = factor(t3_l$subj, levels=c("primeRecon_2", "primeRecon_4", 
                                       "primeRecon_5", "primeRecon_6",
                                       "primeRecon_7", "primeRecon_8", 
                                       "primeRecon_9", "primeRecon_10",
                                       "primeRecon_11", "primeRecon_13",
                                       "primeRecon_14", "primeRecon_15",
                                       "primeRecon_16", "primeRecon_17",
                                       "primeRecon_18", "primeRecon_19",
                                       "primeRecon_20", "primeRecon_21",
                                       "primeRecon_22", "primeRecon_23",
                                       "primeRecon_24"), 
                   labels=c(2,4,5,6,7,8,9,10,11,13,14,15,16,17,18,19,20,21,22,23,24))
t3_l<-t3_l %>% filter(subj!=99)

t3_l$rseCon = factor(t3_l$rseCon, levels=c("init","rep")) 
t3_l$Roi = factor(t3_l$hRoi, levels=c("r_ppa","r_ffa"), labels=c("ppa","ffa"))
t3_l <- t3_l %>% select(subj, Roi, rseCon, z)

# lower ROIs
t4 <- read.csv("data/rse_summary_ms_nsm.csv", header = T)
glimpse(t4, width = 70)
## Rows: 300
## Columns: 9
## $ subj   <chr> "primeRecon_2", "primeRecon_2", "primeRecon_2", "prim…
## $ bgCon  <chr> "rFFA", "rFFA", "rFFA", "rFFA", "rFFA", "rFFA", "rFFA…
## $ fcCon  <chr> "low", "low", "low", "low", "low", "high", "high", "h…
## $ rseCon <chr> "init", "rep", "rse", "n_rs", "n_re", "init", "rep", …
## $ v1     <dbl> 1.07381400, 1.20381300, -0.12999900, 212.00000000, 28…
## $ v2     <dbl> 0.8155636000, 0.7989575000, 0.0166061900, 347.0000000…
## $ v3     <dbl> 0.84924320, 0.82558880, 0.02365439, 318.00000000, 306…
## $ v4     <dbl> 1.23035000, 1.14499500, 0.08535463, 172.00000000, 117…
## $ v_all  <dbl> 0.930998800, 0.939097800, -0.008098876, 949.000000000…

t4_l <- gather(t4, vRoi, z, v1:v_all, factor_key=TRUE)
t4_l$subj = factor(t4_l$subj, levels=c("primeRecon_2", "primeRecon_4", 
                                       "primeRecon_5", "primeRecon_6",
                                       "primeRecon_7", "primeRecon_8", 
                                       "primeRecon_9", "primeRecon_10",
                                       "primeRecon_11", "primeRecon_13",
                                       "primeRecon_14", "primeRecon_15",
                                       "primeRecon_16", "primeRecon_17",
                                       "primeRecon_18", "primeRecon_19",
                                       "primeRecon_20", "primeRecon_21",
                                       "primeRecon_22", "primeRecon_23",
                                       "primeRecon_24"), 
                   labels=c(2,4,5,6,7,8,9,10,11,13,14,15,16,17,18,19,20,21,22,23,24))
t4_l <- t4_l %>% filter(subj!=99)
t4_l$Roi = factor(t4_l$vRoi) 
t4_l$bgCon = factor(t4_l$bgCon)

t4_l$fcCon =factor(t4_l$fcCon, levels=c("low","high","all"))
t4_l$rseCon =factor(t4_l$rseCon, levels=c("init","rep","rse","n_rs","n_re"))

t4_l <- t4_l %>% select(subj, bgCon, Roi, fcCon, rseCon, z) %>% 
  filter(fcCon=="all", rseCon!="rse", rseCon!="n_rs", rseCon!="n_re")
t4_l <- t4_l %>% select(subj, Roi, rseCon, z)

t34 <- rbind(t3_l, t4_l)
t34$rseCon <- factor(t34$rseCon)

glimpse(t34, width = 70)
## Rows: 280
## Columns: 4
## $ subj   <fct> 2, 2, 2, 2, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7,…
## $ Roi    <fct> ppa, ppa, ffa, ffa, ppa, ppa, ffa, ffa, ppa, ppa, ffa…
## $ rseCon <fct> init, rep, init, rep, init, rep, init, rep, init, rep…
## $ z      <dbl> 0.40584820, 0.18225210, 1.31945300, 1.28147900, 0.224…


각 관심 영역에서 제시 조건에 따른 활성화 값의 요약치는 아래와 같다.


t34_allL.1 <- t34 %>% 
  group_by(subj, Roi, rseCon) %>%
  dplyr::summarise(z=mean(z)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'Roi'. You can override using the `.groups` argument.
# t34_allL.1 %>% kable(digits=2)

# subject-level, wide format
t34_allW.1 <- t34_allL.1 %>% spread(key=rseCon, value = z)

# t34_allW.1 <- t34_allL.1 %>% spread(key=Roi, value = z)
# t34_allW.1 %>% kable(digits=2)

# summary table: grand mean
t34_allG.1 <- t34_allL.1 %>% group_by(Roi, rseCon) %>%
  summarise(z.m = mean(z), z.sd = sd(z)) %>%
  ungroup()
## `summarise()` has grouped output by 'Roi'. You can override using the `.groups` argument.
t34_allG.1$z.se <- Rmisc::summarySEwithin(data = t34_allL.1, measurevar = "z", 
                                         idvar = "subj", withinvars = c("Roi", "rseCon"))$se
t34_allG.1$z.ci <- Rmisc::summarySEwithin(data = t34_allL.1, measurevar = "z", 
                                         idvar = "subj", withinvars = c("Roi", "rseCon"))$ci
t34_allG.1 <- t34_allG.1 %>% 
  mutate(lower.ci = z.m-z.ci,
         upper.ci = z.m+z.ci)
t34_allG.1 %>% select(Roi, rseCon, z.m) %>% spread(key=rseCon, value=z.m) %>%  kable(digits=2)
Roi init rep
ppa 0.39 0.32
ffa 2.34 2.13
v1 2.30 2.10
v2 1.96 1.80
v3 1.96 1.80
v4 2.01 1.86
v_all 2.05 1.88


아래 그래프는 각 관심 영역에서 제시 조건에 따른 활성화 값(z)을 나타낸다. 그래프에서 Initial > Repeated 으로 활성화가 나타나는 것을 Repetition Suppression으로 해석한다. 점은 평균, error bar는 95% ci를 나타내었다.


t34.all.plot1 <- ggplot(data=t34_allL.1, aes(x=Roi, y=z, fill=rseCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_segment(data=filter(t34_allW.1, Roi == "ppa"), inherit.aes = FALSE,
             aes(x=0.8, y=filter(t34_allW.1, Roi == "ppa")$init,
                 xend=1.2, yend=filter(t34_allW.1, Roi == "ppa")$rep),
             color="gray80") +
  geom_segment(data=filter(t34_allW.1, Roi == "ffa"), inherit.aes = FALSE,
               aes(x=1.8, y=filter(t34_allW.1, Roi == "ffa")$init,
                   xend=2.2, yend=filter(t34_allW.1, Roi == "ffa")$rep),
               color="gray80") +
  geom_segment(data=filter(t34_allW.1, Roi == "v1"), inherit.aes = FALSE,
               aes(x=2.8, y=filter(t34_allW.1, Roi == "v1")$init,
                   xend=3.2, yend=filter(t34_allW.1, Roi == "v1")$rep),
             color="gray80") +
  geom_segment(data=filter(t34_allW.1, Roi == "v2"), inherit.aes = FALSE,
               aes(x=3.8, y=filter(t34_allW.1, Roi == "v2")$init,
                   xend=4.2, yend=filter(t34_allW.1, Roi == "v2")$rep),
               color="gray80") +
  geom_segment(data=filter(t34_allW.1, Roi == "v3"), inherit.aes = FALSE,
               aes(x=4.8, y=filter(t34_allW.1, Roi == "v3")$init,
                   xend=5.2, yend=filter(t34_allW.1, Roi == "v3")$rep),
               color="gray80") +
  geom_segment(data=filter(t34_allW.1, Roi == "v4"), inherit.aes = FALSE,
               aes(x=5.8, y=filter(t34_allW.1, Roi == "v4")$init,
                   xend=6.2, yend=filter(t34_allW.1, Roi == "v4")$rep),
               color="gray80") +
  geom_segment(data=filter(t34_allW.1, Roi == "v_all"), inherit.aes = FALSE,
               aes(x=6.8, y=filter(t34_allW.1, Roi == "v_all")$init,
                   xend=7.2, yend=filter(t34_allW.1, Roi == "v_all")$rep),
               color="gray80") +
  geom_point(data=t34_allL.1, aes(x=Roi, y=z, fill=rseCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t34_allG.1, aes(x = Roi, y=z.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("PPA","FFA","V1","V2","V3","V4","All VC")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),  # , "#235796"
                    labels = c("Initial", "Repeated")) +
  coord_cartesian(ylim = c(-1, 5), clip = "on") +
  labs(x = "", y = "Activation (z)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")

#ggsave("rse_main1_plot.jpg", plot = t34.all.plot1, width=14, height=6, unit='in', dpi=600)
t34.all.plot1


각 관심 영역(ffa, ppa, v1, v2, v3, v4, all Vc)에서 자극 제시 조건(Initial vs. Repeated)을 요인으로 6 x 2 혼합 요인 분산분석을 수행하였다. 분석 결과에서 각 관심 영역 간 차이는 관심 효과가 아님에 유의할 필요가 있다.


t34.aov1 <- aov_ez(id="subj", dv="z", data = t34_allL.1, within = c("Roi", "rseCon"))
#summary(t4.aov1)
nice(t34.aov1, es="pes") %>% kable(digits=3)
Effect df MSE F pes p.value
Roi 3.96, 75.28 0.31 80.76 *** .810 <.001
rseCon 1, 19 0.18 9.59 ** .335 .006
Roi:rseCon 2.40, 45.63 0.02 2.88 + .132 .057


각 관심 영역에서 자극 제시 조건에 따른 차이를 사후 검정하였다.


# post-hoc - difference between rgCon, by vROI
t34.aov1.m1 <- emmeans(t34.aov1, pairwise ~ rseCon | Roi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t34.aov1.m1$contrasts %>% kable(digits=3)
contrast Roi estimate SE df t.ratio p.value
init - rep ppa 0.065 0.052 19 1.250 0.226
init - rep ffa 0.208 0.063 19 3.312 0.004
init - rep v1 0.203 0.072 19 2.834 0.011
init - rep v2 0.157 0.055 19 2.825 0.011
init - rep v3 0.152 0.049 19 3.081 0.006
init - rep v4 0.151 0.048 19 3.135 0.005
init - rep v_all 0.166 0.056 19 2.974 0.008
#plot(t4.aov1.m1, horizontal = T, comparison=T)


사후 검정 결과, PPA를 제외하고 모든 영역에서 Repetition Suppression 효과(Initial > Repeated)가 나타났다.


2.3 bgFC 강도에 따른 Repetition Suppression (Median Split)


관심 시각 영역 (V1~V4) 에서 Repetition Suppression Effect가 FFA-VC Voxel-wise Background Connectivity(bgFC)의 강도에 따라 조절되는지 검증하였다. 자극 제시 조건은 Initial vs. Repeated, FFA-VC bgFC 강도는 Median Split을 통해 Low, High로 구분되었다. 이 분석에서는 Voxel별 bgFC 값이 주요한 변인이므로 전처리 과정에서 Spatial Smoothing이 수행되지 않은 데이터를 사용하였다.


t4 <- read.csv("data/rse_summary_ms_nsm.csv", header = T)

t4_l <- gather(t4, vRoi, z, v1:v_all, factor_key=TRUE)
# change class of main factors: double to factor
t4_l$subj = factor(t4_l$subj, levels=c("primeRecon_2", "primeRecon_4", 
                                       "primeRecon_5", "primeRecon_6",
                                       "primeRecon_7", "primeRecon_8", 
                                       "primeRecon_9", "primeRecon_10",
                                       "primeRecon_11", "primeRecon_13",
                                       "primeRecon_14", "primeRecon_15",
                                       "primeRecon_16", "primeRecon_17",
                                       "primeRecon_18", "primeRecon_19",
                                       "primeRecon_20", "primeRecon_21",
                                       "primeRecon_22", "primeRecon_23",
                                       "primeRecon_24"), 
                   labels=c(2,4,5,6,7,8,9,10,11,13,14,15,16,17,18,19,20,21,22,23,24))
t4_l <- t4_l %>% filter(subj!=99)

t4_l$vRoi = factor(t4_l$vRoi) 
t4_l$bgCon = factor(t4_l$bgCon)
t4_l$fcCon =factor(t4_l$fcCon, levels=c("low","high","all"))
t4_l$rseCon =factor(t4_l$rseCon, levels=c("init","rep","rse","n_rs","n_re"))
t4_l <- t4_l %>% select(subj, bgCon, vRoi, fcCon, rseCon, z)

glimpse(t4_l, width = 70)
## Rows: 1,500
## Columns: 6
## $ subj   <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, 4,…
## $ bgCon  <fct> rFFA, rFFA, rFFA, rFFA, rFFA, rFFA, rFFA, rFFA, rFFA,…
## $ vRoi   <fct> v1, v1, v1, v1, v1, v1, v1, v1, v1, v1, v1, v1, v1, v…
## $ fcCon  <fct> low, low, low, low, low, high, high, high, high, high…
## $ rseCon <fct> init, rep, rse, n_rs, n_re, init, rep, rse, n_rs, n_r…
## $ z      <dbl> 1.07381400, 1.20381300, -0.12999900, 212.00000000, 28…


각 관심 영역에서 자극 제시 조건과 bgFC 강도에 따른 활성화 값의 요약치는 아래와 같다.


t4_l_1 <- t4_l %>% filter(fcCon!="all", rseCon!="rse", rseCon!="n_rs", rseCon!="n_re")
t4_allL.5 <- t4_l_1 %>% group_by(subj, bgCon, vRoi, fcCon, rseCon) %>%
  dplyr::summarise(z=mean(z)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'bgCon', 'vRoi', 'fcCon'. You can override using the `.groups` argument.
# t4_allL.5 %>% kable(digits=2)

# subject-level, wide format
t4_allW.5 <- t4_allL.5 %>% spread(key=rseCon, value = z)
# t4_allW.5 %>% kable(digits=2)

# summary table: grand mean
t4_allG.5 <- t4_allL.5 %>% group_by(fcCon, vRoi,  rseCon) %>%
  summarise(z.m = mean(z), z.sd = sd(z)) %>%
  ungroup()
## `summarise()` has grouped output by 'fcCon', 'vRoi'. You can override using the `.groups` argument.
t4_allG.5$z.se <- Rmisc::summarySEwithin(data = t4_allL.5, measurevar = "z", 
                                         idvar = "subj", withinvars = c("fcCon", "vRoi", "rseCon"))$se
t4_allG.5$z.ci <- Rmisc::summarySEwithin(data = t4_allL.5, measurevar = "z", 
                                         idvar = "subj", withinvars = c("fcCon", "vRoi", "rseCon"))$ci
t4_allG.5 <- t4_allG.5 %>% 
  mutate(lower.ci = z.m-z.ci,
         upper.ci = z.m+z.ci)
# for between-subject design. check help(summarySE)
# t4_allG.5 %>% kable(digits=2)
# t4_allG.5 %>% select(fcCon, vRoi, rseCon, z.m) %>% spread(key=rseCon, value=z.m) %>% kable(digits=2)
t4_allG.5 %>% select(vRoi, rseCon, fcCon, z.m) %>% spread(key=rseCon, value=z.m) %>% kable(digits=2)
vRoi fcCon init rep
v1 low 1.23 1.08
v1 high 3.37 3.11
v2 low 1.08 0.97
v2 high 2.83 2.63
v3 low 1.04 0.94
v3 high 2.87 2.67
v4 low 1.10 1.00
v4 high 2.91 2.71
v_all low 1.12 1.00
v_all high 2.98 2.77


아래 그래프 중 A는 각 관심 영역에서 자극 제시 조건과 bgFC 강도에 따른 활성화 값(z), B는 각 관심 영역에서 bgFC 강도와 자극 제시 조건에 따른 활성화 값(z)을 나타낸다. A에서 Initial > Repeated로 활성화가 나타나는 것을 Repetition Suppression으로 해석한다. B는 bgFC 강도에 따른 전반적인 활성화 차이를 나타낸다. 점은 평균, error bar는 95% ci를 나타내었다.


t4.all.plot5 <- ggplot(data=t4_allL.5, aes(x=vRoi, y=z, fill=rseCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  geom_point(data=t4_allL.5, aes(x=vRoi, y=z, fill=rseCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t4_allG.5, aes(x = vRoi, y=z.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~fcCon, scales="free_x", space = "free",
               labeller = labeller(fcCon = c("low"="Low bgFC voxels","high"="High bgFC voxels"))) +
  scale_x_discrete(labels=c("V1","V2","V3","V4","all VC")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),  # , "#235796"
                    labels = c("Initial", "Repeated")) +
  coord_cartesian(ylim = c(-2, 10), clip = "on") +
  labs(x = " ", y = "Activation (z)") +
  # ggtitle("ppa-v connectivity") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.5, 1, 0.5, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# t4.all.plot5
#ggsave("rse_main2_1_ms_plot.jpg", plot = t4.all.plot5, width=10, height=6, unit='in', dpi=600)


t4.all.plot5.1 <- ggplot(data=t4_allL.5, aes(x=rseCon, y=z, fill=fcCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  geom_point(data=t4_allL.5, aes(x=rseCon, y=z, fill=fcCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t4_allG.5, aes(x = rseCon, y=z.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~vRoi, scales="free_x", space = "free", 
             labeller = labeller(vRoi = c("v1" = "V1","v2" = "V2","v3" = "V3","v4" = "V4","v_all" = "all VC"))) +
  scale_x_discrete(labels=c("Initial","Repeated")) +
  scale_fill_manual(values = c("#B1D4E0", "#145DA0"),  # , "#235796" "#B1D4E0", "#2E8BC0","#0C2D48","#145DA0"
                    labels = c("low bgFC voxels", "high bgFC voxels")) +
  coord_cartesian(ylim = c(-2, 10), clip = "on") +
  labs(x = " ", y = "Activation (z)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        # axis.text.x = element_text(face = "plain", vjust=0.2, hjust= 0.5, angle = -35, size = 15, color = "black"), 
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain",size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.5, 1, 0.5, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# t4.all.plot5.1
#ggsave("rse_main2_2_ms_plot.jpg", plot = t4.all.plot5.1, width=10, height=6, unit='in', dpi=600)

ggarrange(t4.all.plot5, t4.all.plot5.1, nrow = 2,
          labels=c("A", "B"))


각 관심 시각 영역(v1, v2, v3, v4)에서 자극 제시 조건(Initial vs. Repeated), bgFC 강도(Low vs. High)를 요인으로 4 x 2 x 2 혼합 요인 분산분석을 수행하였다. 분석 결과에서 각 관심 영역 간 차이는 관심 효과가 아님에 유의할 필요가 있다.


t4.aov5 <- aov_ez(id="subj", dv="z", data=t4_allL.5, within = c("vRoi","fcCon","rseCon"))
#summary(t4.aov5)
nice(t4.aov5, es="pes") %>% kable(digits=3)
Effect df MSE F pes p.value
vRoi 2.31, 43.88 0.46 5.26 ** .217 .007
fcCon 1, 19 1.63 204.78 *** .915 <.001
rseCon 1, 19 0.30 9.22 ** .327 .007
vRoi:fcCon 2.02, 38.44 0.23 3.86 * .169 .029
vRoi:rseCon 1.94, 36.89 0.01 1.59 .077 .217
fcCon:rseCon 1, 19 0.04 7.60 * .286 .013
vRoi:fcCon:rseCon 2.47, 47.02 0.00 0.24 .013 .829


관심 영역, 자극 제시 조건, bgFC 강도의 주효과, 관심 영역과 bgFC 강도 간의 상호작용, bgFC 강도와 자극 제시 조건 간의 상호작용이 유의하였다. 각 관심 영역에서 자극 제시 조건과 bgFC 강도에 따른 차이를 사후 검정하였다.


# post-hoc - difference between rgCon, by bgFC, vROI
t4.aov5.m1 <- emmeans(t4.aov5, pairwise ~ rseCon | fcCon + vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t4.aov5.m1$contrasts %>% kable(digits=3)
contrast fcCon vRoi estimate SE df t.ratio p.value
init - rep low v1 0.144 0.064 19 2.236 0.038
init - rep high v1 0.263 0.085 19 3.080 0.006
init - rep low v2 0.109 0.043 19 2.531 0.020
init - rep high v2 0.204 0.073 19 2.813 0.011
init - rep low v3 0.099 0.036 19 2.741 0.013
init - rep high v3 0.206 0.065 19 3.158 0.005
init - rep low v4 0.102 0.033 19 3.119 0.006
init - rep high v4 0.201 0.066 19 3.053 0.007
init - rep low v_all 0.115 0.043 19 2.687 0.015
init - rep high v_all 0.217 0.072 19 3.012 0.007
# plot(t4.aov5.m1, horizontal = T, comparison=T)


모든 시각 영역에서 high bgFC 조건의 repetition suppression이 유의하게 나타났다.


t4.aov5.m2 <- emmeans(t4.aov5, pairwise ~ fcCon | rseCon + vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t4.aov5.m2$contrasts %>% kable(digits=3)
contrast rseCon vRoi estimate SE df t.ratio p.value
low - high init v1 -2.142 0.134 19 -15.959 0
low - high rep v1 -2.023 0.129 19 -15.667 0
low - high init v2 -1.747 0.144 19 -12.163 0
low - high rep v2 -1.651 0.140 19 -11.802 0
low - high init v3 -1.834 0.158 19 -11.630 0
low - high rep v3 -1.727 0.155 19 -11.121 0
low - high init v4 -1.810 0.167 19 -10.859 0
low - high rep v4 -1.712 0.163 19 -10.501 0
low - high init v_all -1.865 0.134 19 -13.930 0
low - high rep v_all -1.763 0.130 19 -13.557 0
# plot(t4.aov5.m1, horizontal = T, comparison=T)


전반적으로, low bgFC 보다 high bgFC의 활성화 값이 높게 나타났다.


2.4 bgFC 강도에 따른 RSE effect (Median Split)


본 연구의 주요 가설은 Repetition Suppression/Enhancement Effect (RSE effect)의 강도가 bgFC 강도에 따라 조절된다는 것이다. 여기서 RSE effect의 강도는 Initial - Repeated의 Difference Score로 정의하였다. RS가 양수이면 Repetition Suppression, 음수이면 Repetition Enhancement를 나타낸다.


관심 시각 영역과 bgFC 강도에 따른 RS effect의 강도의 요약치는 아래와 같다.



t4_l_2 <- t4_l %>% filter(rseCon!="n_rs", rseCon!="n_re", rseCon!="init", 
                          rseCon!="rep", fcCon!="all")
t4_W <- t4_l_2 %>% spread(key=fcCon, value=z)
# t4_W %>% kable(digits=2)

t4_l.dif <- t4_l_2
# t4_l.dif %>% kable(digits=2)

# summary table: grand mean
t4_difG <- t4_l.dif %>% group_by(bgCon, vRoi, fcCon) %>%
  summarise(z.m = mean(z), z.sd = sd(z)) %>%
  ungroup()
## `summarise()` has grouped output by 'bgCon', 'vRoi'. You can override using the `.groups` argument.
t4_difG$z.se <- Rmisc::summarySEwithin(data = t4_l.dif, measurevar = "z", 
                                       idvar = "subj", withinvars = c("bgCon","fcCon", "vRoi"))$se
t4_difG$z.ci <- Rmisc::summarySEwithin(data = t4_l.dif, measurevar = "z", 
                                       idvar = "subj", withinvars = c("bgCon","fcCon", "vRoi"))$ci
t4_difG <- t4_difG %>% 
  mutate(lower.ci = z.m-z.ci,
         upper.ci = z.m+z.ci)
t4_difG %>% select(vRoi, fcCon, z.m) %>% spread(fcCon, z.m) %>%  kable(digits=2)
vRoi low high
v1 0.14 0.26
v2 0.11 0.20
v3 0.10 0.21
v4 0.10 0.20
v_all 0.12 0.22


아래 그래프에는 RSE effect를 관심 시각 영역과 bgFC 강도에 따라 나타내었다. 점은 평균, error bar는 95% CI를 나타낸다.


t4_l.dif.4.1 <- t4_l.dif
t4_w.dif.4.1 <- t4_l.dif.4.1 %>% spread(key=fcCon, value=z)
# t4_w.dif.4.1 %>% kable(digits=2)
t4_difG.4.1 <- t4_difG
# t4_difG.4.1 %>% kable(digits=2)
t4.diff.plot4.1 <- ggplot(data=t4_l.dif.4.1, aes(x=vRoi, y=z, fill=fcCon)) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  stat_summary(fun = mean, geom = "bar", position="dodge", 
               na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_point(data=t4_l.dif.4.1, aes(x=vRoi, y=z, fill=fcCon), 
             position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_segment(data=filter(t4_w.dif.4.1, vRoi == "v1"), inherit.aes = FALSE,
               aes(x=0.8, y=filter(t4_w.dif.4.1, vRoi == "v1")$low,
                   xend=1.2, yend=filter(t4_w.dif.4.1, vRoi == "v1")$high),
               color="gray80") +
  geom_segment(data=filter(t4_w.dif.4.1, vRoi == "v2"), inherit.aes = FALSE,
               aes(x=1.8, y=filter(t4_w.dif.4.1, vRoi == "v2")$low,
                   xend=2.2, yend=filter(t4_w.dif.4.1, vRoi == "v2")$high),
               color="gray80") +
  geom_segment(data=filter(t4_w.dif.4.1, vRoi == "v3"), inherit.aes = FALSE,
               aes(x=2.8, y=filter(t4_w.dif.4.1, vRoi == "v3")$low,
                   xend=3.2, yend=filter(t4_w.dif.4.1, vRoi == "v3")$high),
               color="gray80") +
  geom_segment(data=filter(t4_w.dif.4.1, vRoi == "v4"), inherit.aes = FALSE,
               aes(x=3.8, y=filter(t4_w.dif.4.1, vRoi == "v4")$low,
                   xend=4.2, yend=filter(t4_w.dif.4.1, vRoi == "v4")$high),
               color="gray80") +
  geom_segment(data=filter(t4_w.dif.4.1, vRoi == "v_all"), inherit.aes = FALSE,
               aes(x=4.8, y=filter(t4_w.dif.4.1, vRoi == "v_all")$low,
                   xend=5.2, yend=filter(t4_w.dif.4.1, vRoi == "v_all")$high),
               color="gray80") +
  geom_pointrange(data=t4_difG.4.1, aes(x = vRoi, y=z.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("V1","V2","V3","V4", "All VC")) +
  scale_fill_manual(values = c("#AFDBF5", "#051422"),  # , "#235796"
                    # "#B1D4E0", "#2E8BC0","#0C2D48","#145DA0"
                    # "#AFDBF5", "#195779", "#5BAED9", "#051422"),  # , "#235796"
                    labels = c("Low bgFC", "High bgFC")) +
  coord_cartesian(ylim = c(-1, 2), clip = "on") +
  labs(x = " ", y = "RSE Effect (z)\n[Initial - Repeated]") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.text.x = element_text(face = "plain", vjust = -2, hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
#ggsave("rse_main3_ms_plot.jpg", plot = t4.diff.plot4.1, width=10, height=6, unit='in', dpi=600)
t4.diff.plot4.1


각 관심 시각 영역(v1, v2, v3, v4, all VC)에서 bgFC 강도(Low vs. High)를 요인으로 5 x 2 혼합 요인 분산분석을 수행하였다. 분석 결과에서 각 관심 영역 간 차이는 관심 효과가 아님에 유의할 필요가 있다.


t4_l.dif.4.1 <-  t4_l.dif
t4.2.aov4.1 <- aov_ez(id="subj", dv="z", data = t4_l.dif.4.1, within = c("vRoi", "fcCon"))
# summary(t4.2.aov4.1)
# anova(t4.2.aov4.1, es="pes")
nice(t4.2.aov4.1, es="pes") %>% kable(digits=2)
Effect df MSE F pes p.value
vRoi 1.94, 36.89 0.02 1.59 .077 .217
fcCon 1, 19 0.07 7.60 * .286 .013
vRoi:fcCon 2.47, 47.02 0.01 0.24 .013 .829


bgFC 강도의 주효과가 유의하였다. 각 관심 영역에서 bgFC 강도에 따른 Difference Score의 차이를 사후 검정하였다.


# post-hoc - same-diff RG - difference between fcCon, by vROI
t4.2.aov4.1.m1 <- emmeans(t4.2.aov4.1, pairwise ~ fcCon | vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t4.2.aov4.1.m1$contrasts %>% kable(digits=3)
contrast vRoi estimate SE df t.ratio p.value
low - high v1 -0.119 0.048 19 -2.504 0.022
low - high v2 -0.096 0.044 19 -2.160 0.044
low - high v3 -0.107 0.036 19 -2.946 0.008
low - high v4 -0.098 0.038 19 -2.577 0.018
low - high v_all -0.102 0.040 19 -2.565 0.019
# plot(t4.2.aov4.1.m1, horizontal = T, comparison=T)


모든 영역에서 high bgFC의 RSE effect가 low bgFC 보다 높게 나타났다. 또한 이 값들은 양수로 나타났다. 이는 모든 영역에서 Repetition Suppression이 나타난다는 것을 보여준다.


2.5 bgFC 강도에 따른 Repetition Suppression (4 Quantile Split)


관심 시각 영역 (V1~V4) 에서 Repetition Suppression Effect가 FFA-VC Voxel-wise Background Connectivity(bgFC)의 강도에 따라 조절되는지 검증하였다. 자극 제시 조건은 Initial vs. Repeated, PPA-VC bgFC 강도는 4 Quantile Split을 통해 Low, Low-med, High-med, High로 구분되었다. 이 분석에서는 Voxel별 bgFC 값이 주요한 변인이므로 전처리 과정에서 Spatial Smoothing이 수행되지 않은 데이터를 사용하였다.


t6 <- read.csv("data/rse_summary_qs_nsm.csv", header = T)
t6_l <- gather(t6, vRoi, z, v1:v_all, factor_key=TRUE)

# change class of main factors: double to factor
t6_l$subj = factor(t6_l$subj, levels=c("primeRecon_2", "primeRecon_4", 
                                       "primeRecon_5", "primeRecon_6",
                                       "primeRecon_7", "primeRecon_8", 
                                       "primeRecon_9", "primeRecon_10",
                                       "primeRecon_11", "primeRecon_13",
                                       "primeRecon_14", "primeRecon_15",
                                       "primeRecon_16", "primeRecon_17",
                                       "primeRecon_18", "primeRecon_19",
                                       "primeRecon_20", "primeRecon_21",
                                       "primeRecon_22", "primeRecon_23",
                                       "primeRecon_24"), 
                   labels=c(2,4,5,6,7,8,9,10,11,13,14,15,16,17,18,19,20,21,22,23,24))
t6_l <- t6_l %>% filter(subj!=99)

t6_l$vRoi = factor(t6_l$vRoi) 
t6_l$bgCon = factor(t6_l$bgCon)
t6_l$fcCon =factor(t6_l$fcCon, levels=c("low","low_med","high_med", "high","all"))
t6_l$rseCon =factor(t6_l$rseCon, levels=c("init","rep","rse","n_rs","n_re"))
t6_l <- t6_l %>% select(subj, bgCon, vRoi, fcCon, rseCon, z)

glimpse(t6_l, width = 70)
## Rows: 2,500
## Columns: 6
## $ subj   <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
## $ bgCon  <fct> rFFA, rFFA, rFFA, rFFA, rFFA, rFFA, rFFA, rFFA, rFFA,…
## $ vRoi   <fct> v1, v1, v1, v1, v1, v1, v1, v1, v1, v1, v1, v1, v1, v…
## $ fcCon  <fct> low, low, low, low, low, low_med, low_med, low_med, l…
## $ rseCon <fct> init, rep, rse, n_rs, n_re, init, rep, rse, n_rs, n_r…
## $ z      <dbl> 0.80972600, 1.02135600, -0.21162960, 89.00000000, 162…


각 관심 영역에서 자극 제시 조건과 bgFC 강도에 따른 활성화 값의 요약치는 아래와 같다.


t6_l_1 <- t6_l %>% filter(fcCon!="all", rseCon!="rse", rseCon!="n_rs", rseCon!="n_re")
t6_allL.5 <- t6_l_1 %>%  group_by(subj, bgCon, vRoi, fcCon, rseCon) %>%
  dplyr::summarise(z=mean(z)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'bgCon', 'vRoi', 'fcCon'. You can override using the `.groups` argument.
# t6_allL.5 %>% kable(digits=2)

# subject-level, wide format
t6_allW.5 <- t6_allL.5 %>% spread(key=rseCon, value = z)
# t6_allW.5 %>% kable(digits=2)

# summary table: grand mean
t6_allG.5 <- t6_allL.5 %>% group_by(fcCon, vRoi,  rseCon) %>%
  summarise(z.m = mean(z), z.sd = sd(z)) %>%
  ungroup()
## `summarise()` has grouped output by 'fcCon', 'vRoi'. You can override using the `.groups` argument.
t6_allG.5$z.se <- Rmisc::summarySEwithin(data = t6_allL.5, measurevar = "z", 
                                         idvar = "subj", withinvars = c("fcCon", "vRoi", "rseCon"))$se
t6_allG.5$z.ci <- Rmisc::summarySEwithin(data = t6_allL.5, measurevar = "z", 
                                         idvar = "subj", withinvars = c("fcCon", "vRoi", "rseCon"))$ci
t6_allG.5 <- t6_allG.5 %>% 
  mutate(lower.ci = z.m-z.ci,
         upper.ci = z.m+z.ci)
t6_allG.5 %>% select(vRoi, fcCon, rseCon, z.m) %>% spread(key=rseCon, value=z.m) %>% kable(digits=2)
vRoi fcCon init rep
v1 low 0.76 0.65
v1 low_med 1.70 1.52
v1 high_med 2.63 2.40
v1 high 4.11 3.81
v2 low 0.73 0.63
v2 low_med 1.44 1.32
v2 high_med 2.16 1.99
v2 high 3.50 3.26
v3 low 0.66 0.58
v3 low_med 1.42 1.30
v3 high_med 2.14 1.98
v3 high 3.61 3.36
v4 low 0.75 0.69
v4 low_med 1.45 1.31
v4 high_med 2.17 2.01
v4 high 3.65 3.42
v_all low 0.73 0.64
v_all low_med 1.50 1.37
v_all high_med 2.26 2.08
v_all high 3.70 3.45


아래 그래프 중 A는 각 관심 영역에서 자극 제시 조건과 bgFC 강도에 따른 활성화 값(z), B는 각 관심 영역에서 bgFC 강도와 자극 제시 조건에 따른 활성화 값(z)을 나타낸다. 점은 평균, error bar는 95% ci를 나타내었다.


t6.all.plot5 <- ggplot(data=t6_allL.5, aes(x=vRoi, y=z, fill=rseCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  geom_point(data=t6_allL.5, aes(x=vRoi, y=z, fill=rseCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t6_allG.5, aes(x = vRoi, y=z.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~fcCon, scales="free_x", space = "free") +
  scale_x_discrete(labels=c("V1","V2","V3","V4","all VC")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),  # , "#235796"
                    labels = c("Initial", "Repeated")) +
  coord_cartesian(ylim = c(-2, 10), clip = "on") +
  labs(x = " ", y = "Activation (z)") +
  # ggtitle("ppa-v connectivity") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.5, 1, 0.5, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# t6.all.plot5
#ggsave("rse_main2_1_qs_plot.jpg", plot = t6.all.plot5, width=10, height=6, unit='in', dpi=600)


t6.all.plot5.1 <- ggplot(data=t6_allL.5, aes(x=rseCon, y=z, fill=fcCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  geom_point(data=t6_allL.5, aes(x=rseCon, y=z, fill=fcCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t6_allG.5, aes(x = rseCon, y=z.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~vRoi, scales="free_x", space = "free", 
             labeller = labeller(vRoi = c("v1" = "V1","v2" = "V2","v3" = "V3","v4" = "V4","v_all" = "all VC"))) +
  scale_x_discrete(labels=c("Initial","Repeated")) +
  scale_fill_manual(values = c("#B1D4E0", "#2E8BC0","#0C2D48","#145DA0"),  # , "#235796"
                    labels = c("Lowest", "Low", "High", "Hgihest"))+
  coord_cartesian(ylim = c(-2, 10), clip = "on") +
  labs(x = " ", y = "Activation (z)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        # axis.text.x = element_text(face = "plain", vjust=0.2, hjust= 0.5, angle = -35, size = 15, color = "black"), 
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain",size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.5, 1, 0.5, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# t6.all.plot5.1
#ggsave("rse_main2_2_qs_plot.jpg", plot = t6.all.plot5.1, width=10, height=6, unit='in', dpi=600)

ggarrange(t6.all.plot5, t6.all.plot5.1, nrow = 2,
          labels=c("A", "B"))


각 관심 시각 영역(v1, v2, v3, v4)에서 자극 제시 조건(Initial vs. Repeated), bgFC 강도(Low vs. Low-Med vs. High-Med vs. High)를 요인으로 4 x 2 x 4 혼합 요인 분산분석을 수행하였다. 분석 결과에서 각 관심 영역 간 차이는 관심 효과가 아님에 유의할 필요가 있다.


t6.aov5 <- aov_ez(id="subj", dv="z", data=t6_allL.5, within = c("vRoi","fcCon","rseCon"))
# summary(t6.aov5)
nice(t6.aov5, es="pes") %>% kable(digits=2)
Effect df MSE F pes p.value
vRoi 2.31, 43.89 0.92 5.26 ** .217 .007
fcCon 1.12, 21.22 4.53 182.10 *** .906 <.001
rseCon 1, 19 0.60 9.22 ** .327 .007
vRoi:fcCon 2.99, 56.87 0.49 3.03 * .137 .037
vRoi:rseCon 1.94, 36.90 0.02 1.60 .077 .217
fcCon:rseCon 1.10, 20.90 0.09 7.65 * .287 .010
vRoi:fcCon:rseCon 4.33, 82.18 0.01 0.64 .032 .649


자극 제시 조건, bgFC 강도의 주효과, 관심 영역과 bgFC 강도 간의 상호작용, 관심 영역, 자극 제시 조건, bgFC 강도의 삼원 상호 작용이 유의하였다. 각 관심 영역에서 자극 제시 조건과 bgFC 강도에 따른 차이를 사후 검정하였다.


t6.aov5.m1 <- emmeans(t6.aov5, pairwise ~ rseCon | fcCon + vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t6.aov5.m1$contrasts %>% kable(digits=3)
contrast fcCon vRoi estimate SE df t.ratio p.value
init - rep low v1 0.113 0.061 19 1.853 0.079
init - rep low_med v1 0.175 0.071 19 2.469 0.023
init - rep high_med v1 0.228 0.078 19 2.938 0.008
init - rep high v1 0.298 0.094 19 3.154 0.005
init - rep low v2 0.097 0.040 19 2.455 0.024
init - rep low_med v2 0.120 0.049 19 2.460 0.024
init - rep high_med v2 0.163 0.062 19 2.606 0.017
init - rep high v2 0.246 0.085 19 2.910 0.009
init - rep low v3 0.083 0.032 19 2.612 0.017
init - rep low_med v3 0.115 0.042 19 2.709 0.014
init - rep high_med v3 0.159 0.054 19 2.930 0.009
init - rep high v3 0.253 0.077 19 3.270 0.004
init - rep low v4 0.065 0.032 19 2.010 0.059
init - rep low_med v4 0.139 0.038 19 3.679 0.002
init - rep high_med v4 0.167 0.055 19 3.023 0.007
init - rep high v4 0.234 0.077 19 3.031 0.007
init - rep low v_all 0.094 0.039 19 2.423 0.026
init - rep low_med v_all 0.136 0.049 19 2.774 0.012
init - rep high_med v_all 0.181 0.062 19 2.892 0.009
init - rep high v_all 0.253 0.082 19 3.070 0.006
# plot(t6.aov5.m1, horizontal = T, comparison=T)


high bgFC 조건에서 Repetition Suppression이 유의하게 나타났다.


t6.aov5.m2 <- emmeans(t6.aov5, pairwise ~ fcCon | rseCon + vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t6.aov5.m2$contrasts %>% kable(digits=3)
contrast rseCon vRoi estimate SE df t.ratio p.value
low - low_med init v1 -0.933 0.073 19 -12.811 0
low - high_med init v1 -1.867 0.127 19 -14.704 0
low - high init v1 -3.349 0.204 19 -16.412 0
low_med - high_med init v1 -0.934 0.075 19 -12.526 0
low_med - high init v1 -2.416 0.166 19 -14.530 0
high_med - high init v1 -1.482 0.117 19 -12.694 0
low - low_med rep v1 -0.872 0.059 19 -14.764 0
low - high_med rep v1 -1.753 0.115 19 -15.188 0
low - high rep v1 -3.165 0.192 19 -16.471 0
low_med - high_med rep v1 -0.881 0.075 19 -11.730 0
low_med - high rep v1 -2.293 0.165 19 -13.885 0
high_med - high rep v1 -1.412 0.114 19 -12.398 0
low - low_med init v2 -0.710 0.070 19 -10.134 0
low - high_med init v2 -1.430 0.117 19 -12.245 0
low - high init v2 -2.775 0.230 19 -12.048 0
low_med - high_med init v2 -0.719 0.063 19 -11.508 0
low_med - high init v2 -2.065 0.193 19 -10.678 0
high_med - high init v2 -1.345 0.143 19 -9.388 0
low - low_med rep v2 -0.687 0.064 19 -10.752 0
low - high_med rep v2 -1.364 0.111 19 -12.267 0
low - high rep v2 -2.626 0.223 19 -11.780 0
low_med - high_med rep v2 -0.677 0.062 19 -10.954 0
low_med - high rep v2 -1.939 0.188 19 -10.312 0
high_med - high rep v2 -1.262 0.137 19 -9.196 0
low - low_med init v3 -0.756 0.058 19 -13.118 0
low - high_med init v3 -1.476 0.118 19 -12.470 0
low - high init v3 -2.949 0.249 19 -11.852 0
low_med - high_med init v3 -0.720 0.073 19 -9.804 0
low_med - high init v3 -2.193 0.212 19 -10.321 0
high_med - high init v3 -1.473 0.154 19 -9.583 0
low - low_med rep v3 -0.724 0.060 19 -12.151 0
low - high_med rep v3 -1.399 0.119 19 -11.777 0
low - high rep v3 -2.779 0.244 19 -11.404 0
low_med - high_med rep v3 -0.675 0.073 19 -9.206 0
low_med - high rep v3 -2.055 0.210 19 -9.789 0
high_med - high rep v3 -1.380 0.151 19 -9.113 0
low - low_med init v4 -0.701 0.076 19 -9.188 0
low - high_med init v4 -1.423 0.134 19 -10.592 0
low - high init v4 -2.899 0.262 19 -11.061 0
low_med - high_med init v4 -0.722 0.076 19 -9.472 0
low_med - high init v4 -2.198 0.226 19 -9.741 0
high_med - high init v4 -1.476 0.162 19 -9.110 0
low - low_med rep v4 -0.627 0.072 19 -8.741 0
low - high_med rep v4 -1.321 0.130 19 -10.148 0
low - high rep v4 -2.731 0.259 19 -10.551 0
low_med - high_med rep v4 -0.695 0.072 19 -9.591 0
low_med - high rep v4 -2.104 0.217 19 -9.681 0
high_med - high rep v4 -1.409 0.157 19 -8.972 0
low - low_med init v_all -0.775 0.052 19 -14.922 0
low - high_med init v_all -1.530 0.098 19 -15.601 0
low - high init v_all -2.975 0.213 19 -13.958 0
low_med - high_med init v_all -0.755 0.057 19 -13.267 0
low_med - high init v_all -2.200 0.185 19 -11.902 0
high_med - high init v_all -1.445 0.134 19 -10.752 0
low - low_med rep v_all -0.733 0.048 19 -15.146 0
low - high_med rep v_all -1.443 0.094 19 -15.409 0
low - high rep v_all -2.816 0.208 19 -13.565 0
low_med - high_med rep v_all -0.710 0.055 19 -12.941 0
low_med - high rep v_all -2.083 0.180 19 -11.560 0
high_med - high rep v_all -1.373 0.132 19 -10.417 0
# plot(t6.aov5.m2, horizontal = T, comparison=T)


전반적으로, low bgFC 에서 high bgFC로 갈수록 활성화 값이 높게 나타났다.


2.6 bgFC 강도에 따른 RSE effect (4 Quantile Split)


Median Split에서와 동일한 분석을 수행하였다.


관심 시각 영역과 bgFC 강도에 따른 RSE effect 강도의 요약치는 아래와 같다.


t6_l_2 <- t6_l %>% filter(rseCon!="n_rs", rseCon!="n_re", rseCon!="init", 
                          rseCon!="rep", fcCon!="all")

t6_W <- t6_l_2 %>% spread(key=fcCon, value=z)
# t6_W %>% kable(digits=2)
t6_l.dif <- t6_l_2
# t6_l.dif %>% kable(digits=2)

# summary table: grand mean
t6_difG <- t6_l.dif %>% group_by(bgCon, vRoi, fcCon) %>%
  summarise(z.m = mean(z), z.sd = sd(z)) %>%
  ungroup()
## `summarise()` has grouped output by 'bgCon', 'vRoi'. You can override using the `.groups` argument.
t6_difG$z.se <- Rmisc::summarySEwithin(data = t6_l.dif, measurevar = "z", 
                                       idvar = "subj", withinvars = c("bgCon","fcCon", "vRoi"))$se
t6_difG$z.ci <- Rmisc::summarySEwithin(data = t6_l.dif, measurevar = "z", 
                                       idvar = "subj", withinvars = c("bgCon","fcCon", "vRoi"))$ci
t6_difG <- t6_difG %>% 
  mutate(lower.ci = z.m-z.ci,
         upper.ci = z.m+z.ci)
t6_difG %>% select(vRoi, fcCon, z.m) %>% spread(fcCon, value=z.m)%>% kable(digits=2)
vRoi low low_med high_med high
v1 0.11 0.17 0.23 0.30
v2 0.10 0.12 0.16 0.25
v3 0.08 0.12 0.16 0.25
v4 0.07 0.14 0.17 0.23
v_all 0.09 0.14 0.18 0.25


아래 그래프에는 RSE effect를 관심 시각 영역과 bgFC 강도에 따라 나타내었다. 점은 평균, error bar는 95% CI를 나타낸다.


t5_l.dif.4.1 <- t6_l.dif 
t5_w.dif.4.1 <- t5_l.dif.4.1 %>% spread(key=fcCon, value=z)
t5_difG.4.1 <- t6_difG
# t5_difG.4.1 %>% kable(digits=2)
t5.diff.plot5.1 <- ggplot(data=t5_l.dif.4.1, aes(x=vRoi, y=z, fill=fcCon)) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  stat_summary(fun = mean, geom = "bar", position="dodge", 
               na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_point(data=t5_l.dif.4.1, aes(x=vRoi, y=z, fill=fcCon), 
             position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v1"), inherit.aes = FALSE,
               aes(x=0.7, y=filter(t5_w.dif.4.1, vRoi == "v1")$low,
                   xend=0.9, yend=filter(t5_w.dif.4.1, vRoi == "v1")$low_med),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v1"), inherit.aes = FALSE,
               aes(x=0.9, y=filter(t5_w.dif.4.1, vRoi == "v1")$low_med,
                   xend=1.1, yend=filter(t5_w.dif.4.1, vRoi == "v1")$high_med),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v1"), inherit.aes = FALSE,
               aes(x=1.1, y=filter(t5_w.dif.4.1, vRoi == "v1")$high_med,
                   xend=1.3, yend=filter(t5_w.dif.4.1, vRoi == "v1")$high),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v2"), inherit.aes = FALSE,
               aes(x=1.7, y=filter(t5_w.dif.4.1, vRoi == "v2")$low,
                   xend=1.9, yend=filter(t5_w.dif.4.1, vRoi == "v2")$low_med),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v2"), inherit.aes = FALSE,
               aes(x=1.9, y=filter(t5_w.dif.4.1, vRoi == "v2")$low_med,
                   xend=2.1, yend=filter(t5_w.dif.4.1, vRoi == "v2")$high_med),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v2"), inherit.aes = FALSE,
               aes(x=2.1, y=filter(t5_w.dif.4.1, vRoi == "v2")$high_med,
                   xend=2.3, yend=filter(t5_w.dif.4.1, vRoi == "v2")$high),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v3"), inherit.aes = FALSE,
               aes(x=2.7, y=filter(t5_w.dif.4.1, vRoi == "v3")$low,
                   xend=2.9, yend=filter(t5_w.dif.4.1, vRoi == "v3")$low_med),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v3"), inherit.aes = FALSE,
               aes(x=2.9, y=filter(t5_w.dif.4.1, vRoi == "v3")$low_med,
                   xend=3.1, yend=filter(t5_w.dif.4.1, vRoi == "v3")$high_med),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v3"), inherit.aes = FALSE,
               aes(x=3.1, y=filter(t5_w.dif.4.1, vRoi == "v3")$high_med,
                   xend=3.3, yend=filter(t5_w.dif.4.1, vRoi == "v3")$high),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v4"), inherit.aes = FALSE,
               aes(x=3.7, y=filter(t5_w.dif.4.1, vRoi == "v4")$low,
                   xend=3.9, yend=filter(t5_w.dif.4.1, vRoi == "v4")$low_med),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v4"), inherit.aes = FALSE,
               aes(x=3.9, y=filter(t5_w.dif.4.1, vRoi == "v4")$low_med,
                   xend=4.1, yend=filter(t5_w.dif.4.1, vRoi == "v4")$high_med),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v4"), inherit.aes = FALSE,
               aes(x=4.1, y=filter(t5_w.dif.4.1, vRoi == "v4")$high_med,
                   xend=4.3, yend=filter(t5_w.dif.4.1, vRoi == "v4")$high),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v_all"), inherit.aes = FALSE,
             aes(x=4.7, y=filter(t5_w.dif.4.1, vRoi == "v_all")$low,
                 xend=4.9, yend=filter(t5_w.dif.4.1, vRoi == "v_all")$low_med),
             color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v_all"), inherit.aes = FALSE,
               aes(x=4.9, y=filter(t5_w.dif.4.1, vRoi == "v_all")$low_med,
                   xend=5.1, yend=filter(t5_w.dif.4.1, vRoi == "v_all")$high_med),
               color="gray80") +
  geom_segment(data=filter(t5_w.dif.4.1, vRoi == "v_all"), inherit.aes = FALSE,
               aes(x=5.1, y=filter(t5_w.dif.4.1, vRoi == "v_all")$high_med,
                   xend=5.3, yend=filter(t5_w.dif.4.1, vRoi == "v_all")$high),
               color="gray80") +
  geom_pointrange(data=t5_difG.4.1, aes(x = vRoi, y=z.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("V1","V2","V3","V4", "All VC")) +
  scale_fill_manual(values = c("#AFDBF5", "#195779", "#5BAED9", "#051422"),  # , "#235796"
                    labels = c("Lowest", "Low",  "High",  "Highest")) +
  coord_cartesian(ylim = c(-1, 2), clip = "on") +
  labs(x = " ", y = "RSE effect (z)\n[Same - Different]") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.text.x = element_text(face = "plain", vjust=-2, hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
        # legend.position=c(0.8, 0.85))
t5.diff.plot5.1

#ggsave("rse_main3_3_ qs_plot.jpg", plot = t5.diff.plot5.1, width=8, height=6, unit='in', dpi=600)


각 관심 시각 영역(v1, v2, v3, v4, all VC)에서 bgFC 강도(Low vs. Low-med vs. High-Med vs. High)를 요인으로 5 x 4 혼합 요인 분산분석을 수행하였다. 분석 결과에서 각 관심 영역 간 차이는 관심 효과가 아님에 유의할 필요가 있다.


t5.2.aov4.1 <- aov_ez(id="subj", dv="z", data = t5_l.dif.4.1, within = c("vRoi", "fcCon"))
# summary(t4.2.aov4.1)
nice(t5.2.aov4.1, es="pes") %>% kable(digits=3)
Effect df MSE F pes p.value
vRoi 1.94, 36.90 0.05 1.60 .077 .217
fcCon 1.10, 20.90 0.18 7.65 * .287 .010
vRoi:fcCon 4.33, 82.18 0.01 0.64 .032 .649


관심 영역의 주효과가 유의하였다. 관심 영역과 bgFC 강도 간의 상호 작용은 유의하지 않았으나, 높은 트렌드를 보였다. 각 관심 영역에서 bgFC 강도에 따른 Difference Score의 차이를 사후 검정하였다.


t5.2.aov4.1.m1 <- emmeans(t5.2.aov4.1, pairwise ~ fcCon | vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t5.2.aov4.1.m1$contrasts %>% kable(digits=3)
contrast vRoi estimate SE df t.ratio p.value
low - low_med v1 -0.061 0.030 19 -2.028 0.213
low - high_med v1 -0.114 0.053 19 -2.166 0.169
low - high v1 -0.184 0.071 19 -2.590 0.078
low_med - high_med v1 -0.053 0.026 19 -2.082 0.195
low_med - high v1 -0.123 0.046 19 -2.665 0.067
high_med - high v1 -0.070 0.027 19 -2.597 0.077
low - low_med v2 -0.023 0.023 19 -0.985 0.759
low - high_med v2 -0.065 0.044 19 -1.488 0.464
low - high v2 -0.149 0.067 19 -2.225 0.152
low_med - high_med v2 -0.043 0.024 19 -1.793 0.307
low_med - high v2 -0.126 0.049 19 -2.567 0.081
high_med - high v2 -0.083 0.031 19 -2.655 0.068
low - low_med v3 -0.032 0.020 19 -1.603 0.401
low - high_med v3 -0.077 0.035 19 -2.224 0.153
low - high v3 -0.170 0.058 19 -2.951 0.038
low_med - high_med v3 -0.044 0.017 19 -2.655 0.068
low_med - high v3 -0.138 0.042 19 -3.296 0.018
high_med - high v3 -0.093 0.029 19 -3.241 0.020
low - low_med v4 -0.074 0.026 19 -2.849 0.046
low - high_med v4 -0.102 0.039 19 -2.594 0.077
low - high v4 -0.169 0.057 19 -2.967 0.036
low_med - high_med v4 -0.028 0.023 19 -1.198 0.635
low_med - high v4 -0.094 0.044 19 -2.143 0.175
high_med - high v4 -0.067 0.027 19 -2.428 0.105
low - low_med v_all -0.042 0.022 19 -1.883 0.268
low - high_med v_all -0.086 0.040 19 -2.138 0.177
low - high v_all -0.159 0.061 19 -2.607 0.075
low_med - high_med v_all -0.045 0.020 19 -2.286 0.137
low_med - high v_all -0.117 0.042 19 -2.819 0.049
high_med - high v_all -0.072 0.026 19 -2.806 0.051
# plot(t4.2.aov4.1.m1, horizontal = T, comparison=T)


v3, v4에서 high bgFC의 Difference Score가 low bgFC 보다 높게 나타났다.


2.7 각 복셀별 RSE effect 와 bgFC 간의 상관 관계


관심 시각 영역 (V1~V4, all VC) 에서 RSE Effect가 FFA-VC Voxel-wise Background Connectivity(bgFC)의 강도에 따라 조절되는지 검증하기 위해 더 직접적인 상관 분석을 수행하였다. 각 참가자 내의 관심 영역에서 voxel들을 전체 voxel, RS voxel, RE voxel의 세 조건으로 나누었다. 각 조건에서 voxel 별로 Initial - Repeated의 Difference Score를 계산한 후 각 voxel 별 FFA-VC bgFC 값과 상관을 계산하였다. RSE는 initial - Repeated로 계산되었기 때문에, 양수는 Repetition Suppression, 음수 값음 Repetition Enhancement 를 의미한다. 전체 복셀에서 0 보다 유의하게 상관계수가 높다는 것은 bgFC 강도가 높을수록 RSE effect의 강도도 높다는 것을 나타낸다. RS 복셀에서 0보다 유의하게 상관계수가 높다는 것은 bgFC 강도가 클수록 Repetition Suppression 효과가 강하다는 것을 나타낸다. 반대로, RE 복셀에서는 음수를 양수로 변환하였다. 따라서 RE 복셀에서 0보다 높은 상관계수는 bgFC 강도가 클수록 Repetition Enhancement 효과가 강하다는 것을 나타낸다.


s1 <- read.csv("data/rse_summary_corr_nsm.csv", header = T)

s1_l <- s1
s1_l$subj = factor(s1_l$subj, levels=c("primeRecon_2", "primeRecon_4", 
                                       "primeRecon_5", "primeRecon_6",
                                       "primeRecon_7", "primeRecon_8", 
                                       "primeRecon_9", "primeRecon_10",
                                       "primeRecon_11", "primeRecon_13",
                                       "primeRecon_14", "primeRecon_15",
                                       "primeRecon_16", "primeRecon_17",
                                       "primeRecon_18", "primeRecon_19",
                                       "primeRecon_20", "primeRecon_21",
                                       "primeRecon_22", "primeRecon_23",
                                       "primeRecon_24"), 
                   labels=c(2,4,5,6,7,8,9,10,11,13,14,15,16,17,18,19,20,21,22,23,24))
s1_l <- s1_l %>% filter(subj!=99)
s1_l$vRoi = factor(s1_l$vRoi, levels=c("v1","v2","v3","v4","v_all")) 
s1_l$rseCon <- factor(s1_l$rseCon, levels=c("all", "rs", "re"), labels=c("all", "RS", "RE"))
glimpse(s1_l, width = 70)
## Rows: 300
## Columns: 5
## $ subj   <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, 4,…
## $ vRoi   <fct> v1, v2, v3, v4, v_all, v1, v2, v3, v4, v_all, v1, v2,…
## $ Zr     <dbl> 0.484906900, 0.035958040, 0.189455700, 0.186666700, 0…
## $ pval   <dbl> 0.000000e+00, 1.921995e-01, 2.574792e-11, 8.000843e-0…
## $ rseCon <fct> all, all, all, all, all, RS, RS, RS, RS, RS, RE, RE, …


각 관심 영역에서 voxel 별 difference score와 bgFC의 상관 관계를 Fisher의 Z transform하여 요약하였다.


s1_allL <- s1_l %>% 
  select(subj, vRoi, rseCon, Zr) %>% 
  group_by(subj, vRoi, rseCon) %>%
  dplyr::summarise(Zr=mean(Zr)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'vRoi'. You can override using the `.groups` argument.
# s1_allL %>% kable(digits=2)

# subject-level, wide format
s1_allW <- s1_allL %>% select(subj, vRoi, rseCon, Zr) %>% spread(key=vRoi, value = Zr)
#s1_allW %>% kable(digits=2)

# summary table: grand mean
s1_allG <- s1_allL %>% group_by(vRoi, rseCon) %>%
  summarise(Zr.m = mean(Zr), Zr.sd = sd(Zr)) %>%
  ungroup()
## `summarise()` has grouped output by 'vRoi'. You can override using the `.groups` argument.
s1_allG$Zr.se <- Rmisc::summarySEwithin(data = s1_allL, measurevar = "Zr", 
                                        idvar = "subj", withinvars = c("vRoi","rseCon"))$se
s1_allG$Zr.ci <- Rmisc::summarySEwithin(data = s1_allL, measurevar = "Zr", 
                                        idvar = "subj", withinvars = c("vRoi","rseCon"))$ci
s1_allG <- s1_allG %>% 
  mutate(lower.ci = Zr.m-Zr.ci,
         upper.ci = Zr.m+Zr.ci)
s1_allG %>% select(rseCon, vRoi, Zr.m) %>% spread(vRoi, Zr.m) %>% kable(digits=2)
rseCon v1 v2 v3 v4 v_all
all 0.15 0.14 0.17 0.15 0.14
RS 0.11 0.10 0.15 0.13 0.12
RE -0.13 -0.07 -0.07 -0.04 -0.07


추가로 개별 참가자의 상관분석에서 획득된 p-value를 나타내었다. p-value는 1 - p-value로 계산되었으므로, 0.95 보다 높은 값은 해당 상관 관계의 p-value가 유의하게 나타났다는 것을 보여준다.


s1_allL1 <- s1_l %>% mutate(pval = 1-pval)

s1_allL1 <- s1_allL1 %>% 
  select(subj, vRoi, rseCon, pval) %>% 
  group_by(subj, vRoi, rseCon) %>%
  dplyr::summarise(pval=mean(pval)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'vRoi'. You can override using the `.groups` argument.
# s1_allL1 %>% kable(digits=2)

s1_allW1 <- s1_allL1 %>% select(subj, vRoi, rseCon, pval) %>% spread(key=vRoi, value = pval)
# s1_allW1 %>% kable(digits=2)

s1_allG1 <- s1_allL1 %>% group_by(vRoi, rseCon) %>%
  summarise(pval.m = mean(pval), pval.sd = sd(pval)) %>%
  ungroup()
## `summarise()` has grouped output by 'vRoi'. You can override using the `.groups` argument.
s1_allG1$pval.se <- Rmisc::summarySEwithin(data = s1_allL1, measurevar = "pval", 
                                         idvar = "subj", withinvars = c("vRoi","rseCon"))$se
s1_allG1$pval.ci <- Rmisc::summarySEwithin(data = s1_allL1, measurevar = "pval", 
                                         idvar = "subj", withinvars = c("vRoi","rseCon"))$ci
s1_allG1 <- s1_allG1 %>% 
  mutate(lower.ci = pval.m-pval.ci,
         upper.ci = pval.m+pval.ci)
s1_allG1 %>% select(rseCon, vRoi, pval.m) %>% spread(vRoi, pval.m) %>% kable(digits=2)
rseCon v1 v2 v3 v4 v_all
all 0.99 0.88 0.89 0.81 0.94
RS 0.88 0.79 0.88 0.91 0.91
RE 0.74 0.76 0.82 0.44 0.85


각 관심 영역의 전체 voxel 조건에서 voxel 별 difference score와 bgFC의 Zr과 p-value 들을 그래프로 나타내었다. 이 분석에는 모든 RSE 복셀이 포함되었다. A는 voxel 별 difference score와 bgFC의 Zr, B는 voxel 별 difference score와 bgFC correlation의 p-value를 나타낸다. 점은 평군, error bar는 95% CI를 나타낸다.


s1_allL.1 <- s1_allL %>% filter(rseCon=="all")
s1_allG.1 <- s1_allG %>% filter(rseCon=="all")
s1.all.plot1 <- ggplot(data=s1_allL.1, aes(x=vRoi, y=Zr, fill=vRoi)) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02, show.legend=F) +
  geom_point(data=s1_allL.1, aes(x=vRoi, y=Zr, fill=vRoi), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=s1_allG.1, aes(x = vRoi, y=Zr.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("V1","V2","V3","V4", "All VC")) +
  scale_fill_manual(values = c("#B1D4E0", "#2E8BC0","#0C2D48","#145DA0","#DBDBDB")) +
  coord_cartesian(ylim = c(-0.1, 0.5), clip = "on") +
  labs(x = "", y = "Correlation (Zr)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# s1.all.plot1
#ggsave("rg_corr_nsm_plot1.jpg", plot = s1.all.plot1, width=8, height=8, unit='in', dpi=600)

s1_allL1.1 <- s1_allL1 %>% filter(rseCon=="all")
s1_allG1.1 <- s1_allG1 %>% filter(rseCon=="all")
s1.all.plot2 <- ggplot(data=s1_allL1.1, aes(x=vRoi, y=pval, fill=vRoi)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02, show.legend=F) +
  geom_hline(yintercept=0.95, linetype='dashed', color='black', alpha =1, size=1) +
  geom_point(data=s1_allL1.1, aes(x=vRoi, y=pval, fill=vRoi), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=s1_allG1.1, aes(x = vRoi, y=pval.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("V1","V2","V3","V4", "All VC")) +
  scale_fill_manual(values = c("#B1D4E0", "#2E8BC0","#0C2D48","#145DA0","#DBDBDB")) +
  coord_cartesian(ylim = c(0.4, 1.1), clip = "on") +
  labs(x = "", y = "1-pvalue of Correlation (RG & bgFC)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# s1.all.plot2
#ggsave("rg_corr_nsm_plot2.jpg", plot = s1.all.plot2, width=8, height=8, unit='in', dpi=600)

ggarrange(s1.all.plot1, s1.all.plot2, ncol=2,
          labels=c("A", "B"))


각 관심 영역에서 voxel 별 difference score와 bgFC의 Zr이 0보다 유의하게 큰 지 확인하기 위해 단일 표본 t-검정을 수행하였다.


s1.v1 <- s1_allL.1 %>% filter(vRoi=="v1")
s1.v2 <- s1_allL.1 %>% filter(vRoi=="v2")
s1.v3 <- s1_allL.1 %>% filter(vRoi=="v3")
s1.v4 <- s1_allL.1 %>% filter(vRoi=="v4")
s1.vall <- s1_allL.1 %>% filter(vRoi=="v_all")


s1.ttest.v1 <- t.test(Zr ~ 1, data=s1.v1,
                       alternative = c("two.sided"),
                       conf.level=.95)
s1.ttest.v1
## 
##  One Sample t-test
## 
## data:  Zr
## t = 2.5114, df = 19, p-value = 0.02122
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.02558223 0.28155901
## sample estimates:
## mean of x 
## 0.1535706

s1.ttest.v2 <- t.test(Zr ~ 1, data=s1.v2,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.v2
## 
##  One Sample t-test
## 
## data:  Zr
## t = 2.305, df = 19, p-value = 0.03262
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.01252152 0.25982233
## sample estimates:
## mean of x 
## 0.1361719

s1.ttest.v3 <- t.test(Zr ~ 1, data=s1.v3,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.v3
## 
##  One Sample t-test
## 
## data:  Zr
## t = 3.1675, df = 19, p-value = 0.005071
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.05732964 0.28068688
## sample estimates:
## mean of x 
## 0.1690083

s1.ttest.v4 <- t.test(Zr ~ 1, data=s1.v4,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.v4
## 
##  One Sample t-test
## 
## data:  Zr
## t = 2.8957, df = 19, p-value = 0.009266
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.04278593 0.26591728
## sample estimates:
## mean of x 
## 0.1543516

s1.ttest.vall <- t.test(Zr ~ 1, data=s1.vall,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.vall
## 
##  One Sample t-test
## 
## data:  Zr
## t = 2.7134, df = 19, p-value = 0.01379
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.03272208 0.25352614
## sample estimates:
## mean of x 
## 0.1431241


아래의 Results 2 분석에서는 상위 영역인 FFA의 RS voxel과 하위 영역 V3&4의 RE voxel 간의 연결성이 FFA RE - V3&4 RE 간의 연결성보다 높게 나타났다. 이는 Error unit이 상위 영역의 Representation Unit으로 Feedforward signal을 보낸다는 PC의 설명과 맞는 것으로 볼 수 있다.


그러나 Results 1의 분석에서 상위 영역과 높은 연결성을 보일수록 VC의 RE 효과가 적다는 것도 확인하였다. 이 두 결과는 얼핏보면 서로 맞지 않는 것으로 보인다. 이는 Results 1의 분석에서 FFA를 RS/RE로 나누어 살펴보지 않았기 때문일 수 있다. 또 다른 가능한 설명은 FFA - VC 연결성이 RS/RE라는 두 유형의 복셀들에 서로 다르게 작용하는데 복셀 수의 차이 때문에 이 효과가 가려진다는 것이다. RS 복셀들에서 RS의 효과가 클수록 연결성이 강하고, RE 복셀들에서 RE의 효과가 클수록 연결성이 강하다고 하더라도, RS 복셀의 수가 RE 복셀들보다 많기 때문에 Results 1에서의 결과가 나타날 가능성이 있다. Results 1의 해당 분석에서 연결성 값이 드라마틱하게 강하지는 않은 이유가 이렇게 혼재된 관계의 voxel들이 함께 분석되어서일 수 있다. 만약 그렇다면, 복셀들을 사전에 RS/RE로 구분하고 상위 영역과의 연결성과의 관계를 살펴보았을 때 RS의 강도 - 연결성 강도 간의 정적 상관과 RE의 강도 - 연결성 강도의 정적 상관이 모두 관찰될 수 있을 것이다. 아래의 분석에서는 Visual Cortex에서 RS/RE를 구분한 후, 각 복셀들에서 RS effect와 RE effect의 강도와 연결성 간의 관계를 살펴보았다.


동일하게 RS voxel과 RE voxel을 나누어 그래프로 나타내었다. A는 voxel 별 difference score와 bgFC의 Zr, B는 voxel 별 difference score와 bgFC correlation의 p-value를 나타낸다. 점은 평군, error bar는 95% CI를 나타낸다.


s1_allL.2 <- s1_allL %>% filter(rseCon!="all")
s1_allG.2 <- s1_allG %>% filter(rseCon!="all")
s1.all.plot3 <- ggplot(data=s1_allL.2, aes(x=vRoi, y=Zr, fill=rseCon)) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02, show.legend=T) +
  geom_point(data=s1_allL.2, aes(x=vRoi, y=Zr, fill=rseCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=s1_allG.2, aes(x = vRoi, y=Zr.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("V1","V2","V3","V4", "All VC")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),  # , "#235796"
                    labels = c("RSE of RSs", "RSE of REs")) +
  coord_cartesian(ylim = c(-0.6, 0.6), clip = "on") +
  labs(x = "", y = "Correlation (Zr)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
s1.all.plot3

#ggsave("rg_corr_nsm_plot3.jpg", plot = s1.all.plot1, width=8, height=8, unit='in', dpi=600)

s1_allL1.2 <- s1_allL1 %>% filter(rseCon!="all")
s1_allG1.2 <- s1_allG1 %>% filter(rseCon!="all")
s1.all.plot4 <- ggplot(data=s1_allL1.2, aes(x=vRoi, y=pval, fill=rseCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02, show.legend=T) +
  geom_hline(yintercept=0.95, linetype='dashed', color='black', alpha =1, size=1) +
  geom_point(data=s1_allL1.2, aes(x=vRoi, y=pval, fill=rseCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=s1_allG1.2, aes(x = vRoi, y=pval.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("V1","V2","V3","V4", "All VC")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),  # , "#235796"
                    labels = c("RSE of RSs", "RSE of REs")) +
  coord_cartesian(ylim = c(0.4, 1.1), clip = "on") +
  labs(x = "", y = "1-pvalue of Correlation (RG & bgFC)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
s1.all.plot4

#ggsave("rg_corr_nsm_plot4.jpg", plot = s1.all.plot4, width=8, height=8, unit='in', dpi=600)

# ggarrange(s1.all.plot3, s1.all.plot4, ncol=2,
          # labels=c("A", "B"))


RS/RE 간에서 연결성의 유형은 유의하게 달랐다. RS 복셀에서는 RS의 크기가 클수록 bgFC 강도도 컸다. 그러나 RE 복셀에서는 RE의 크기가 적을수록 bgFC의 강도가 컸다. 해석하자면, RS 효과가 클수록 상위 영역과 강한 연결성을 보이고, RE 효과가 클수록 상위 영역과 약한 연결성을 보인다는 것이다. 이는 앞선 가설과 일치하지 않는다. RE가 Error Unit이라고 본다면, 상위 영역으로의 feedforward 연결이 존재할 수 있다. 특히 상위 영역의 RS 복셀들로의 Feedforward 연결이 나타날 수 있다. 그러나 현재 분석에서 RE 효과 자체와 상위 영역과 연결성은 역상관을 보이는 것 같다. 물론 현재 상위 영역인 FFA의 timecourse는 RS/RE를 구분하지 않은 연결성이다. 아래에는 FFA에서 RS, RE를 구분하여 연결성을 계산하여 추가할 예정이다.


s1.aov1.1 <- aov_ez(id="subj", dv="Zr", data = s1_allL.2 , within = c("vRoi", "rseCon"))
nice(s1.aov1.1, es="pes") %>% kable(digits=2)
Effect df MSE F pes p.value
vRoi 2.54, 48.31 0.01 2.24 .105 .105
rseCon 1, 19 0.17 11.16 ** .370 .003
vRoi:rseCon 2.27, 43.05 0.02 1.03 .051 .374


RS/RE 간에서 연결성의 유형은 통계적으로 유의하게 달랐다.


# post-hoc - same-diff RG - difference between fcCon, by vROI
s1.aov1.1.m1 <- emmeans(s1.aov1.1, pairwise ~ rseCon | vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
s1.aov1.1.m1$contrasts %>% kable(digits=3)
contrast vRoi estimate SE df t.ratio p.value
RS - RE v1 0.235 0.067 19 3.501 0.002
RS - RE v2 0.166 0.065 19 2.566 0.019
RS - RE v3 0.217 0.061 19 3.552 0.002
RS - RE v4 0.174 0.067 19 2.584 0.018
RS - RE v_all 0.183 0.059 19 3.085 0.006
# plot(t4.2.aov4.1.m1, horizontal = T, comparison=T)


마찬가지로 각 조건에서 연결성의 강도가 0보다 유의하게 큰지 검정하였다.


먼저 RS voxel들에서의 결과는 다음과 같다.


s1.v1 <- s1_allL.2 %>% filter(rseCon=="RS", vRoi=="v1")
s1.v2 <- s1_allL.2 %>% filter(rseCon=="RS", vRoi=="v2")
s1.v3 <- s1_allL.2 %>% filter(rseCon=="RS", vRoi=="v3")
s1.v4 <- s1_allL.2 %>% filter(rseCon=="RS", vRoi=="v4")
s1.vall <- s1_allL.2 %>% filter(rseCon=="RS", vRoi=="v_all")

s1.ttest.v1 <- t.test(Zr ~ 1, data=s1.v1,
                       alternative = c("two.sided"),
                       conf.level=.95)
s1.ttest.v1
## 
##  One Sample t-test
## 
## data:  Zr
## t = 2.5384, df = 19, p-value = 0.02005
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.01872324 0.19471259
## sample estimates:
## mean of x 
## 0.1067179

s1.ttest.v2 <- t.test(Zr ~ 1, data=s1.v2,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.v2
## 
##  One Sample t-test
## 
## data:  Zr
## t = 2.1949, df = 19, p-value = 0.0408
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.004521744 0.190263047
## sample estimates:
## mean of x 
## 0.0973924

s1.ttest.v3 <- t.test(Zr ~ 1, data=s1.v3,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.v3
## 
##  One Sample t-test
## 
## data:  Zr
## t = 3.0831, df = 19, p-value = 0.006121
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.04704495 0.24595036
## sample estimates:
## mean of x 
## 0.1464977

s1.ttest.v4 <- t.test(Zr ~ 1, data=s1.v4,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.v4
## 
##  One Sample t-test
## 
## data:  Zr
## t = 2.8867, df = 19, p-value = 0.009451
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.03691985 0.23164598
## sample estimates:
## mean of x 
## 0.1342829

s1.ttest.vall <- t.test(Zr ~ 1, data=s1.vall,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.vall
## 
##  One Sample t-test
## 
## data:  Zr
## t = 3.0365, df = 19, p-value = 0.00679
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.03602116 0.19584675
## sample estimates:
## mean of x 
##  0.115934


RE voxel들에서의 결과는 다음과 같다.


s1.v1 <- s1_allL.2 %>% filter(rseCon=="RE", vRoi=="v1")
s1.v2 <- s1_allL.2 %>% filter(rseCon=="RE", vRoi=="v2")
s1.v3 <- s1_allL.2 %>% filter(rseCon=="RE", vRoi=="v3")
s1.v4 <- s1_allL.2 %>% filter(rseCon=="RE", vRoi=="v4")
s1.vall <- s1_allL.2 %>% filter(rseCon=="RE", vRoi=="v_all")

s1.ttest.v1 <- t.test(Zr ~ 1, data=s1.v1,
                       alternative = c("two.sided"),
                       conf.level=.95)
s1.ttest.v1
## 
##  One Sample t-test
## 
## data:  Zr
## t = -3.4884, df = 19, p-value = 0.002459
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.20573451 -0.05143334
## sample estimates:
##  mean of x 
## -0.1285839

s1.ttest.v2 <- t.test(Zr ~ 1, data=s1.v2,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.v2
## 
##  One Sample t-test
## 
## data:  Zr
## t = -2.1775, df = 19, p-value = 0.04225
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.135135825 -0.002672145
## sample estimates:
##   mean of x 
## -0.06890399

s1.ttest.v3 <- t.test(Zr ~ 1, data=s1.v3,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.v3
## 
##  One Sample t-test
## 
## data:  Zr
## t = -2.3604, df = 19, p-value = 0.0291
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.132550355 -0.007956934
## sample estimates:
##   mean of x 
## -0.07025364

s1.ttest.v4 <- t.test(Zr ~ 1, data=s1.v4,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.v4
## 
##  One Sample t-test
## 
## data:  Zr
## t = -1.1974, df = 19, p-value = 0.2459
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.10809801  0.02942383
## sample estimates:
##   mean of x 
## -0.03933709

s1.ttest.vall <- t.test(Zr ~ 1, data=s1.vall,
                      alternative = c("two.sided"),
                      conf.level=.95)
s1.ttest.vall
## 
##  One Sample t-test
## 
## data:  Zr
## t = -2.2267, df = 19, p-value = 0.03826
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.129985433 -0.004022209
## sample estimates:
##   mean of x 
## -0.06700382


2.8 관심 영역에서 RS/RE 복셀의 비율


가설은 상위 영역과 연결성이 강한 복셀일 수록 RSE의 RS 효과가 크게 나타날 수 있다는 것이었다. 만약 그렇다면 RS/RE 복셀의 개수 자체에서도 차이가 나타날 수 있을까? 일반적으로 RS와 RE의 비율은 70 대 30 정도로 나누어져 있는 듯하다. 이것을 연결성 강도에 따라 나누어서 봤을 때, 우리의 가설이 맞다면 연결성이 강할수록, 즉 median split에서 high 조건인 경우에 RS 복셀과 RE 복셀 간 비율 차이가 더 커질 것이다.


# Higher ROIs
t8 <- read.csv("data/rse_summary_ffa_nsm.csv", header = T)
t8_l <- t8

# change class of main factors: double to factor
t8_l$subj = factor(t8_l$subj, levels=c("primeRecon_2", "primeRecon_4",
                                       "primeRecon_5", "primeRecon_6",
                                       "primeRecon_7", "primeRecon_8",
                                       "primeRecon_9", "primeRecon_10",
                                       "primeRecon_11", "primeRecon_13",
                                       "primeRecon_14", "primeRecon_15",
                                       "primeRecon_16", "primeRecon_17",
                                       "primeRecon_18", "primeRecon_19",
                                       "primeRecon_20", "primeRecon_21",
                                       "primeRecon_22", "primeRecon_23",
                                       "primeRecon_24"),
                   labels=c(2,4,5,6,7,8,9,10,11,13,14,15,16,17,18,19,20,21,22,23,24))
t8_l <- t8_l %>% filter(subj!=99)
t8_l$hRoi = factor(t8_l$hRoi) 
t8_l$rsCon =factor(t8_l$rsCon)
t8_l <- t8_l %>% filter(rsCon != "rs_init", rsCon != "re_init",
                        rsCon != "rs_rep", rsCon != "re_rep",
                        rsCon != "rs_rse", rsCon != "re_rse")
t8_l$rsCon =factor(t8_l$rsCon, levels=c("n_rs","n_re"))
t8_l$z <- t8_l$value
t8_l$Roi <- t8_l$hRoi
t8_l$fcCon <- "all"
t8_l <- t8_l %>% select(subj, Roi, fcCon, rsCon, z)
t8_l_3 <- t8_l
t8_l_4 <- t8_l_3 %>% spread(key=rsCon, value=z)
t8_l_4 <- t8_l_4 %>% mutate(rs_v = n_rs/(n_rs+n_re), re_v = n_re/(n_rs+n_re) )
t8_l_4 <- t8_l_4 %>% select(subj, Roi, fcCon, rs_v, re_v)
t8_l_4 <- gather(t8_l_4, rsCon, prop, c("rs_v", "re_v"), factor_key=T)
t8_l_4$rseCon <- t8_l_4$rsCon
t8_l_4 <- t8_l_4 %>% select(subj, Roi, fcCon, rseCon, prop)

# lower ROIs
t4_l_3 <- t4_l %>% filter(rseCon!="init", rseCon!="rep", rseCon!="rse")
t4_l_4 <- t4_l_3 %>% spread(key=rseCon, value=z)
t4_l_4 <- t4_l_4 %>% mutate(rs_v = n_rs/(n_rs+n_re), re_v = n_re/(n_rs+n_re) )
t4_l_4 <- t4_l_4 %>% select(subj, bgCon, vRoi, fcCon, rs_v, re_v)
t4_l_4 <- gather(t4_l_4, rseCon, prop, c("rs_v", "re_v"), factor_key=T)
t4_l_4$Roi <- t4_l_4$vRoi
t4_l_4 <- t4_l_4 %>% select(subj, Roi, fcCon, rseCon, prop)

t84_l_4 <- rbind(t4_l_4, t8_l_4)
t84_l_4$subj <- factor(t84_l_4$subj)
t84_l_4$Roi <- factor(t84_l_4$Roi, levels=c("ffa", "v1", "v2", "v3", "v4", "v_all"))


연결성 강도를 low, high, 전체의 세 가지로 나누어서 각 관심 영역의 RS/RE 복셀 비율을 나타내면 다음과 같다.


t84_1_allL.5 <- t84_l_4 %>% group_by(subj, Roi, fcCon, rseCon) %>%
  dplyr::summarise(prop=mean(prop)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'Roi', 'fcCon'. You can override using the `.groups` argument.

# subject-level, wide format
t84_1_allW.5 <- t84_1_allL.5 %>% spread(key=rseCon, value = prop)

# summary table: grand mean
t84_1_allG.5 <- t84_1_allL.5 %>% group_by(fcCon, Roi, rseCon) %>%
  summarise(prop.m = mean(prop), prop.sd = sd(prop)) %>%
  ungroup()
## `summarise()` has grouped output by 'fcCon', 'Roi'. You can override using the `.groups` argument.
t84_1_allG.5$prop.se <- Rmisc::summarySEwithin(data = t84_1_allL.5, measurevar = "prop", 
                                              idvar = "subj", withinvars = c("fcCon", "Roi", "rseCon"))$se
t84_1_allG.5$prop.ci <- Rmisc::summarySEwithin(data = t84_1_allL.5, measurevar = "prop", 
                                              idvar = "subj", withinvars = c("fcCon", "Roi", "rseCon"))$ci
t84_1_allG.5 <- t84_1_allG.5 %>% 
  mutate(lower.ci = prop.m-prop.ci,
         upper.ci = prop.m+prop.ci,
         lower.se = prop.m-prop.se,
         upper.se = prop.m+prop.se)
# for between-subject design. check help(summarySE)
# t4_1_allG.5 %>% kable(digits=2)
t84_1_allG.5 %>% select(fcCon, Roi, rseCon, prop.m) %>% spread(key=rseCon, value=prop.m) %>% kable(digits=2)
fcCon Roi rs_v re_v
low v1 0.63 0.37
low v2 0.60 0.40
low v3 0.59 0.41
low v4 0.61 0.39
low v_all 0.60 0.40
high v1 0.70 0.30
high v2 0.67 0.33
high v3 0.68 0.32
high v4 0.68 0.32
high v_all 0.68 0.32
all ffa 0.66 0.34
all v1 0.66 0.34
all v2 0.63 0.37
all v3 0.64 0.36
all v4 0.64 0.36
all v_all 0.64 0.36


먼저 연결성을 반영하지 않고, 전체 복셀에서 RS/RE의 비율을 살펴보았다. FFA, V1~V4, All VC에 걸쳐서 대략 7대 3, 6대 4의 비율이 유지되는 것으로 보인다.


t84_1_allL.5.1 <- t84_1_allL.5 %>% filter(fcCon=="all")
t84_1_allG.5.1 <- t84_1_allG.5 %>% filter(fcCon=="all")
t84_1.all.plot5 <- ggplot(data=t84_1_allL.5.1, aes(x=Roi, y=prop, fill=rseCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  geom_point(data=t84_1_allL.5.1, aes(x=Roi, y=prop, fill=rseCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t84_1_allG.5.1, aes(x = Roi, y=prop.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("FFA","V1","V2","V3","V4","all VC")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),  # , "#235796"
                    labels = c("RS voxels", "RE voxels")) +
  coord_cartesian(ylim = c(0, 1), clip = "on") +
  labs(x = " ", y = "Proportion of voxels (%)") +
  # ggtitle("ppa-v connectivity") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.5, 1, 0.5, "cm"), 
        legend.title = element_blank(),
        legend.position=c(0.8, 0.85))

#ggsave("rse_voxels1_ms_plot.jpg", plot = t84_1.all.plot5, width=10, height=6, unit='in', dpi=600)
t84_1.all.plot5


영역 간 차이나 영역, RS/RE 조건 간의 상호작용은 나타나지 않았다.


t84_1.aov5 <- aov_ez(id="subj", dv="prop", data=t84_1_allL.5.1, within = c("Roi","rseCon"))
nice(t84_1.aov5, es="pes")
## Anova Table (Type 3 tests)
## 
## Response: prop
##       Effect          df  MSE        F   pes p.value
## 1        Roi 1.70, 32.39 0.00     0.00 <.001   >.999
## 2     rseCon       1, 19 0.34 14.73 **  .437    .001
## 3 Roi:rseCon 2.49, 47.24 0.03     0.54  .028    .625
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
## 
## Sphericity correction method: GG


사후 검정 결과에서도, 영역 간 RS/RE 복셀 비율은 차이가 없는 것처럼 보인다.


t84_1.aov5.m1 <- emmeans(t84_1.aov5, pairwise ~ Roi | rseCon, adjust = "tukey") # adjust="bon" , "tukey", "none"
t84_1.aov5.m1$contrasts %>% kable(digits=4)
contrast rseCon estimate SE df t.ratio p.value
ffa - v1 rs_v 0.0005 0.0373 19 0.0121 1.0000
ffa - v2 rs_v 0.0314 0.0349 19 0.8993 0.9420
ffa - v3 rs_v 0.0270 0.0347 19 0.7786 0.9679
ffa - v4 rs_v 0.0198 0.0385 19 0.5149 0.9949
ffa - v_all rs_v 0.0220 0.0336 19 0.6532 0.9850
v1 - v2 rs_v 0.0309 0.0202 19 1.5343 0.6479
v1 - v3 rs_v 0.0266 0.0246 19 1.0778 0.8842
v1 - v4 rs_v 0.0194 0.0265 19 0.7319 0.9753
v1 - v_all rs_v 0.0215 0.0153 19 1.4030 0.7249
v2 - v3 rs_v -0.0044 0.0152 19 -0.2867 0.9997
v2 - v4 rs_v -0.0116 0.0224 19 -0.5150 0.9949
v2 - v_all rs_v -0.0094 0.0086 19 -1.0927 0.8783
v3 - v4 rs_v -0.0072 0.0202 19 -0.3553 0.9991
v3 - v_all rs_v -0.0050 0.0112 19 -0.4502 0.9973
v4 - v_all rs_v 0.0021 0.0176 19 0.1217 1.0000
ffa - v1 re_v -0.0005 0.0373 19 -0.0121 1.0000
ffa - v2 re_v -0.0314 0.0349 19 -0.8993 0.9420
ffa - v3 re_v -0.0270 0.0347 19 -0.7786 0.9679
ffa - v4 re_v -0.0198 0.0385 19 -0.5149 0.9949
ffa - v_all re_v -0.0220 0.0336 19 -0.6532 0.9850
v1 - v2 re_v -0.0309 0.0202 19 -1.5343 0.6479
v1 - v3 re_v -0.0266 0.0246 19 -1.0778 0.8842
v1 - v4 re_v -0.0194 0.0265 19 -0.7319 0.9753
v1 - v_all re_v -0.0215 0.0153 19 -1.4030 0.7249
v2 - v3 re_v 0.0044 0.0152 19 0.2867 0.9997
v2 - v4 re_v 0.0116 0.0224 19 0.5150 0.9949
v2 - v_all re_v 0.0094 0.0086 19 1.0927 0.8783
v3 - v4 re_v 0.0072 0.0202 19 0.3553 0.9991
v3 - v_all re_v 0.0050 0.0112 19 0.4502 0.9973
v4 - v_all re_v -0.0021 0.0176 19 -0.1217 1.0000
# plot(t84_1.aov5.m1, horipropontal = T, comparison=T)


2.9 bgFC 강도에 따른 관심 영역의 RS/RE 복셀 비율 차이


bgFC 강도를 Median Split하여 Low/High로 구분한 후, 각 관심 시각 영역에서의 RS/RE 복셀 비율이 bgFC 강도와 상호작용하는지 살펴보았다.


t4_l_3 <- t4_l %>% filter(rseCon!="init", rseCon!="rep", rseCon!="rse")
t4_l_4 <- t4_l_3 %>% spread(key=rseCon, value=z)
t4_l_4 <- t4_l_4 %>% mutate(rs_v = n_rs/(n_rs+n_re), re_v = n_re/(n_rs+n_re) )
t4_l_4 <- t4_l_4 %>% select(subj, bgCon, vRoi, fcCon, rs_v, re_v)
t4_l_4 <- gather(t4_l_4, rseCon, prop, c("rs_v", "re_v"), factor_key=T)


요약치는 다음과 같다.


t4_1_allL.5 <- t4_l_4 %>% filter(fcCon!="all") %>% group_by(subj, bgCon, vRoi, fcCon, rseCon) %>%
  dplyr::summarise(prop=mean(prop)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'bgCon', 'vRoi', 'fcCon'. You can override using the `.groups` argument.
# t4_1_allL.5 %>% kable(digits=2)

# subject-level, wide format
t4_1_allW.5 <- t4_1_allL.5 %>% spread(key=rseCon, value = prop)
# t4_1_allW.5 %>% kable(digits=2)

# summary table: grand mean
t4_1_allG.5 <- t4_1_allL.5 %>% group_by(fcCon, vRoi,  rseCon) %>%
  summarise(prop.m = mean(prop), prop.sd = sd(prop)) %>%
  ungroup()
## `summarise()` has grouped output by 'fcCon', 'vRoi'. You can override using the `.groups` argument.
t4_1_allG.5$prop.se <- Rmisc::summarySEwithin(data = t4_1_allL.5, measurevar = "prop", 
                                              idvar = "subj", withinvars = c("fcCon", "vRoi", "rseCon"))$se
t4_1_allG.5$prop.ci <- Rmisc::summarySEwithin(data = t4_1_allL.5, measurevar = "prop", 
                                              idvar = "subj", withinvars = c("fcCon", "vRoi", "rseCon"))$ci
t4_1_allG.5 <- t4_1_allG.5 %>% 
  mutate(lower.ci = prop.m-prop.ci,
         upper.ci = prop.m+prop.ci,
         lower.se = prop.m-prop.se,
         upper.se = prop.m+prop.se)
t4_1_allG.5 %>% select(vRoi, fcCon, rseCon, prop.m) %>% spread(key=rseCon, value=prop.m) %>% kable(digits=2)
vRoi fcCon rs_v re_v
v1 low 0.63 0.37
v1 high 0.70 0.30
v2 low 0.60 0.40
v2 high 0.67 0.33
v3 low 0.59 0.41
v3 high 0.68 0.32
v4 low 0.61 0.39
v4 high 0.68 0.32
v_all low 0.60 0.40
v_all high 0.68 0.32


A는 각 bgFC 강도에서 시각 영역과 RS/RE비율 간 상호작용을 나타내고, B는 각 시각 영역에서 RS/RE 비율과 bgFC 강도 간의 상호작용을 나타낸다. 점은 평균, errorbar는 95% CI를 표시하였다.


t4_1_allL.5.2 <- t4_1_allL.5 %>% filter(fcCon!="all")
t4_1_allG.5.2 <- t4_1_allG.5 %>% filter(fcCon!="all")
t4_1.all.plot5.2 <- ggplot(data=t4_1_allL.5.2, aes(x=vRoi, y=prop, fill=rseCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  geom_point(data=t4_1_allL.5.2, aes(x=vRoi, y=prop, fill=rseCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t4_1_allG.5.2, aes(x = vRoi, y=prop.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~fcCon, scales="free_x", space = "free",
             labeller = labeller(fcCon = c("low"="Low bgFC voxels","high"="High bgFC voxels"))) +
  scale_x_discrete(labels=c("V1","V2","V3","V4","all VC")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),  # , "#235796"
                    labels = c("RS voxels", "RE voxels")) +
  coord_cartesian(ylim = c(0, 1.2), clip = "on") +
  labs(x = " ", y = "Proportion of voxels (%)") +
  # ggtitle("ppa-v connectivity") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.5, 1, 0.5, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# t4_1.all.plot5.2
#ggsave("rse_voxels2_ms_plot.jpg", plot = t4_1.all.plot5.2, width=10, height=6, unit='in', dpi=600)

t4_1_allL.6 <- t4_1_allL.5 %>% filter(fcCon!="all")
t4_1_allG.6 <- t4_1_allG.5 %>% filter(fcCon!="all")
t4_1.all.plot5.1 <- ggplot(data=t4_1_allL.6, aes(x=rseCon, y=prop, fill=fcCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=0, linetype='dashed', color='gray80', alpha =1, size=1) +
  geom_point(data=t4_1_allL.6, aes(x=rseCon, y=prop, fill=fcCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t4_1_allG.6, aes(x = rseCon, y=prop.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~vRoi, scales="free_x", space = "free", 
             labeller = labeller(vRoi = c("v1" = "V1","v2" = "V2","v3" = "V3","v4" = "V4","v_all" = "all VC"))) +
  scale_x_discrete(labels=c("RS","RE")) +
  scale_fill_manual(values = c("#B1D4E0", "#145DA0"),  # , "#235796" "#B1D4E0", "#2E8BC0","#0C2D48","#145DA0"
                    labels = c("low bgFC voxels", "high bgFC voxels")) +
  coord_cartesian(ylim = c(0, 1.2), clip = "on") +
  labs(x = " ", y = "Proportion of voxels (%)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        # axis.text.x = element_text(face = "plain", vjust=0.2, hjust= 0.5, angle = -35, size = 15, color = "black"), 
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain",size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.5, 1, 0.5, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# t4_1.all.plot5.1
#ggsave("rse_voxels3_ms_plot.jpg", plot = t4_1.all.plot5.1, width=10, height=6, unit='in', dpi=600)

ggarrange(t4_1.all.plot5.2, t4_1.all.plot5.1, ncol = 2, labels=c("A", "B"))


시각 관심 영역(V1~V4, all VC), bgFC 강도(low, high), RS/RE 복셀비율 (RS 비율, RE 비율)을 요인으로 5 x 2 x 2 혼합 요인 분산분석을 수행하였다. bgFC 강도와 RS/RE 복셀 비율 간의 상호작용이 유의하였다.


t4_1.aov5 <- aov_ez(id="subj", dv="prop", data=t4_1_allL.6, within = c("vRoi","fcCon","rseCon"))
# summary(t4_1.aov5)
nice(t4_1.aov5, es="pes") %>% kable(digits=2)
Effect df MSE F pes p.value
vRoi 1.57, 29.88 0.00 0.00 <.001 >.999
fcCon 1, 19 -0.00 -0.00 <.001 >.999
rseCon 1, 19 0.59 13.73 ** .419 .002
vRoi:fcCon 1.61, 30.56 -0.00 -0.00 <.001 >.999
vRoi:rseCon 2.43, 46.21 0.02 0.79 .040 .484
fcCon:rseCon 1, 19 0.05 9.52 ** .334 .006
vRoi:fcCon:rseCon 2.39, 45.46 0.01 0.10 .005 .933


가설과 일치하게, 연결성이 강한 조건이 약한 조건보다 RS 비율과 RE 비율 간의 차이가 크게 나타나는 것으로 보인다.


t4_1.aov5.m1 <- emmeans(t4_1.aov5, pairwise ~ rseCon | fcCon + vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t4_1.aov5.m1$contrasts %>% kable(digits=3)
contrast fcCon vRoi estimate SE df t.ratio p.value
rs_v - re_v low v1 0.253 0.097 19 2.618 0.017
rs_v - re_v high v1 0.393 0.107 19 3.673 0.002
rs_v - re_v low v2 0.192 0.069 19 2.799 0.011
rs_v - re_v high v2 0.331 0.096 19 3.438 0.003
rs_v - re_v low v3 0.189 0.062 19 3.061 0.006
rs_v - re_v high v3 0.351 0.088 19 3.976 0.001
rs_v - re_v low v4 0.216 0.062 19 3.463 0.003
rs_v - re_v high v4 0.353 0.091 19 3.867 0.001
rs_v - re_v low v_all 0.210 0.069 19 3.048 0.007
rs_v - re_v high v_all 0.351 0.093 19 3.784 0.001
# plot(t4_1.aov5.m1, horipropontal = T, comparison=T)


t4_1.aov5.m2 <- emmeans(t4_1.aov5, pairwise ~ fcCon | rseCon + vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t4_1.aov5.m2$contrasts %>% kable(digits=3)
contrast rseCon vRoi estimate SE df t.ratio p.value
low - high rs_v v1 -0.070 0.035 19 -2.016 0.058
low - high re_v v1 0.070 0.035 19 2.016 0.058
low - high rs_v v2 -0.070 0.030 19 -2.291 0.034
low - high re_v v2 0.070 0.030 19 2.291 0.034
low - high rs_v v3 -0.081 0.023 19 -3.529 0.002
low - high re_v v3 0.081 0.023 19 3.529 0.002
low - high rs_v v4 -0.069 0.022 19 -3.162 0.005
low - high re_v v4 0.069 0.022 19 3.162 0.005
low - high rs_v v_all -0.070 0.025 19 -2.808 0.011
low - high re_v v_all 0.070 0.025 19 2.808 0.011
# plot(t4_1.aov5.m2, horipropontal = T, comparison=T)


2.10 각 관심 영역 조건과 bgFC 조건에서 bgFC의 평균값


전반적인 연결성 정도를 확인하기 위해 각 관심 영역 조건에서 bgFC 조건의 connectivity 평균값을 표시하였다.


# median split
s2 <- read.csv("data/bg_summary_ms_nsm.csv", header = T)
glimpse(s2, width = 70)
## Rows: 100
## Columns: 5
## $ subj    <chr> "primeRecon_2", "primeRecon_2", "primeRecon_2", "pri…
## $ vRoi    <chr> "v1", "v2", "v3", "v4", "v_all", "v1", "v2", "v3", "…
## $ low_Zr  <dbl> 0.077025535, 0.040950922, 0.037788743, 0.050176767, …
## $ high_Zr <dbl> 0.2376537, 0.1971796, 0.1958559, 0.2048766, 0.210664…
## $ all_Zr  <dbl> 0.15641723, 0.11834508, 0.11616207, 0.12677080, 0.12…
# check number of trials for each condition/subj
table(s2$subj)
## 
## primeRecon_10 primeRecon_11 primeRecon_13 primeRecon_14 primeRecon_15 
##             5             5             5             5             5 
## primeRecon_16 primeRecon_17 primeRecon_18 primeRecon_19  primeRecon_2 
##             5             5             5             5             5 
## primeRecon_21 primeRecon_22 primeRecon_23 primeRecon_24  primeRecon_4 
##             5             5             5             5             5 
##  primeRecon_5  primeRecon_6  primeRecon_7  primeRecon_8  primeRecon_9 
##             5             5             5             5             5
table(s2$vRoi, s2$subj) 
##        
##         primeRecon_10 primeRecon_11 primeRecon_13 primeRecon_14 primeRecon_15
##   v_all             1             1             1             1             1
##   v1                1             1             1             1             1
##   v2                1             1             1             1             1
##   v3                1             1             1             1             1
##   v4                1             1             1             1             1
##        
##         primeRecon_16 primeRecon_17 primeRecon_18 primeRecon_19 primeRecon_2
##   v_all             1             1             1             1            1
##   v1                1             1             1             1            1
##   v2                1             1             1             1            1
##   v3                1             1             1             1            1
##   v4                1             1             1             1            1
##        
##         primeRecon_21 primeRecon_22 primeRecon_23 primeRecon_24 primeRecon_4
##   v_all             1             1             1             1            1
##   v1                1             1             1             1            1
##   v2                1             1             1             1            1
##   v3                1             1             1             1            1
##   v4                1             1             1             1            1
##        
##         primeRecon_5 primeRecon_6 primeRecon_7 primeRecon_8 primeRecon_9
##   v_all            1            1            1            1            1
##   v1               1            1            1            1            1
##   v2               1            1            1            1            1
##   v3               1            1            1            1            1
##   v4               1            1            1            1            1

s2_l <- gather(s2, fcCon, fc, low_Zr:all_Zr, factor_key=T)


# change class of main factors: double to factor
s2_l$subj = factor(s2_l$subj, levels=c("primeRecon_2", "primeRecon_4", 
                                       "primeRecon_5", "primeRecon_6",
                                       "primeRecon_7", "primeRecon_8", 
                                       "primeRecon_9", "primeRecon_10",
                                       "primeRecon_11", "primeRecon_13",
                                       "primeRecon_14", "primeRecon_15",
                                       "primeRecon_16", "primeRecon_17",
                                       "primeRecon_18", "primeRecon_19",
                                       "primeRecon_20", "primeRecon_21",
                                       "primeRecon_22", "primeRecon_23",
                                       "primeRecon_24"), 
                   labels=c(2,4,5,6,7,8,9,10,11,13,14,15,16,17,18,19,20,21,22,23,24))
s2_l <- s2_l %>% filter(subj!=99)
s2_l$vRoi = factor(s2_l$vRoi, levels=c("v1","v2","v3","v4","v_all")) 
s2_l$fcCon = factor(s2_l$fcCon, levels=c("low_Zr", "high_Zr", "all_Zr"), 
                    labels=c("low","high","all")) 

# 4 quantile
s3 <- read.csv("data/bg_summary_qs_nsm.csv", header = T)
s3_l <- gather(s3, fcCon, fc, low_Zr:all_Zr, factor_key=T)
# change class of main factors: double to factor
s3_l$subj = factor(s3_l$subj, levels=c("primeRecon_2", "primeRecon_4", 
                                       "primeRecon_5", "primeRecon_6",
                                       "primeRecon_7", "primeRecon_8", 
                                       "primeRecon_9", "primeRecon_10",
                                       "primeRecon_11", "primeRecon_13",
                                       "primeRecon_14", "primeRecon_15",
                                       "primeRecon_16", "primeRecon_17",
                                       "primeRecon_18", "primeRecon_19",
                                       "primeRecon_20", "primeRecon_21",
                                       "primeRecon_22", "primeRecon_23",
                                       "primeRecon_24"), 
                   labels=c(2,4,5,6,7,8,9,10,11,13,14,15,16,17,18,19,20,21,22,23,24))
s3_l <- s3_l %>% filter(subj!=99)
s3_l$vRoi = factor(s3_l$vRoi, levels=c("v1","v2","v3","v4","v_all")) 
s3_l$fcCon = factor(s3_l$fcCon, levels=c("low_Zr", "low_med_Zr","high_med_Zr","high_Zr", "all_Zr"), 
                    labels=c("low","low_med","high_med","high","all")) 


Median split 과 4 Quantile 버전에서 각각의 요약치는 아래와 같다


Median Split 버전


# subject-level, long format
s2_allL <- s2_l %>% filter(fcCon!="all") %>% 
  group_by(subj, vRoi, fcCon) %>%
  dplyr::summarise(fc=mean(fc)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'vRoi'. You can override using the `.groups` argument.
# s2_allL %>% kable(digits=2)

# subject-level, wide format
s2_allW <- s2_allL %>% spread(key=fcCon, value = fc)
# s2_allW %>% kable(digits=2)

# summary table: grand mean
s2_allG <- s2_allL %>% group_by(vRoi, fcCon) %>%
  summarise(fc.m = mean(fc), fc.sd = sd(fc)) %>%
  ungroup()
## `summarise()` has grouped output by 'vRoi'. You can override using the `.groups` argument.
s2_allG$fc.se <- Rmisc::summarySEwithin(data = s2_allL, measurevar = "fc", 
                                        idvar = "subj", withinvars = c("vRoi", "fcCon"))$se
s2_allG$fc.ci <- Rmisc::summarySEwithin(data = s2_allL, measurevar = "fc", 
                                        idvar = "subj", withinvars = c("vRoi", "fcCon"))$ci
s2_allG <- s2_allG %>% 
  mutate(lower.ci = fc.m-fc.ci,
         upper.ci = fc.m+fc.ci)
s2_allG %>% select(vRoi, fcCon, fc.m) %>% spread(key=fcCon, value=fc.m) %>%  kable(digits=2)
vRoi low high
v1 0.05 0.20
v2 0.04 0.20
v3 0.05 0.20
v4 0.06 0.21
v_all 0.05 0.20


4 Quantile 버전

# subject-level, long format
s3_allL <- s3_l %>% filter(fcCon!="all") %>% 
  group_by(subj, vRoi, fcCon) %>%
  dplyr::summarise(fc=mean(fc)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'vRoi'. You can override using the `.groups` argument.
# s3_allL %>% kable(digits=2)

# subject-level, wide format
s3_allW <- s3_allL %>% spread(key=fcCon, value = fc)
# s3_allW %>% kable(digits=2)

# summary table: grand mean
s3_allG <- s3_allL %>% group_by(vRoi, fcCon) %>%
  summarise(fc.m = mean(fc), fc.sd = sd(fc)) %>%
  ungroup()
## `summarise()` has grouped output by 'vRoi'. You can override using the `.groups` argument.
s3_allG$fc.se <- Rmisc::summarySEwithin(data = s3_allL, measurevar = "fc", 
                                        idvar = "subj", withinvars = c("vRoi", "fcCon"))$se
s3_allG$fc.ci <- Rmisc::summarySEwithin(data = s3_allL, measurevar = "fc", 
                                        idvar = "subj", withinvars = c("vRoi", "fcCon"))$ci
s3_allG <- s3_allG %>% 
  mutate(lower.ci = fc.m-fc.ci,
         upper.ci = fc.m+fc.ci)
s3_allG %>% select(vRoi, fcCon, fc.m) %>% spread(key=fcCon, value=fc.m) %>%  kable(digits=2)
vRoi low low_med high_med high
v1 0.01 0.09 0.15 0.24
v2 0.00 0.08 0.15 0.25
v3 0.01 0.08 0.15 0.26
v4 0.03 0.10 0.16 0.27
v_all 0.01 0.09 0.15 0.26


그래포로 나타나면 아래와 같다.


s2.all.plot1 <- ggplot(data=s2_allL, aes(x=vRoi, y=fc, fill=fcCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_point(data=s2_allL, aes(x=vRoi, y=fc, fill=fcCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=s2_allG, aes(x = vRoi, y=fc.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("V1","V2","V3","V4","all VC")) +
  scale_fill_manual(values = c("#AFDBF5", "#051422"),  # , "#235796"
                    labels = c("Low",  "High")) +
  
  coord_cartesian(ylim = c(-0.2, 0.5), clip = "on") +
  labs(x = "", y = "BgFC (Zr), by Median Split") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# s2.all.plot1
#ggsave("bg_ms_plot1.jpg", plot = s2.all.plot1, width=8, height=8, unit='in', dpi=600)

s3.all.plot1 <- ggplot(data=s3_allL, aes(x=vRoi, y=fc, fill=fcCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_point(data=s3_allL, aes(x=vRoi, y=fc, fill=fcCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=s3_allG, aes(x = vRoi, y=fc.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("V1","V2","V3","V4","all VC")) +
  scale_fill_manual(values = c("#AFDBF5", "#195779", "#5BAED9", "#051422"),  # , "#235796"
                    labels = c("Lowest", "Low",  "High",  "Highest")) +
  coord_cartesian(ylim = c(-0.2, 0.5), clip = "on") +
  labs(x = "", y = "BgFC (Zr), by 4 Quantile Split") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# s3.all.plot1
#ggsave("bg_qs_plot2.jpg", plot = s3.all.plot1, width=8, height=8, unit='in', dpi=600)

s23.plot1 <- ggarrange(s2.all.plot1, s3.all.plot1, nrow = 2)
#ggsave("bg_qsms_plot.jpg", plot = s23.plot1, width=10, height=8, unit='in', dpi=600)
s23.plot1


3 Results 2 - Visual RS/RE voxel - FFA RS/RE voxel 간의 bgFC


Predictive Coding Model을 RSE에 적용한다고 했을 때, RS는 Representation unit, RE는 Error unit에 대응될 수 있다. PC에서 Error 유닛은 상위 영역의 Representation unit에 Feedforward signal을 전송한다. 단순하게 연결성의 맥락으로 옮겨서보자면, 하위 영역의 RE와 상위 영역의 RS 간 연결성이 강하게 나타날 수 있다. (그렇다면 상위 영역에서의 prediction, feedback signal은 어디에서 오는가..?)


상위 영역인 FFA를 RS/RE로 구분하고 하위 영역인 VC의 RS/RE와의 연결성을 분석해볼 수 있다. 이 경우 2 x 2의 4가지 연결성이 계산될 수 있다. 만약 PC를 그대로 적용한 가설이 맞는다면, VC의 RE와 FFA의 RS간 연결성이 다른 연결성보다 강하게 나타날 것이다. VC RE - FFA RE, VC RS - FFA RS, VC RS - FFA RE는 어떻게 될 것인가? 현재 이 세가지 비교에 대해서는 특별한 가설이 없다. 일단 확인해보자.


+ 추가로 해볼 것, 이전 분석(Results 1)에서는 전체 FFA의 time Course를 활용했다. 이것을 RS/RE로 나눠본다면 어떨까?


3.1 Visual RS/RE voxel과 상위 영역 RS/RE voxel 간의 연결성


상위 영역인 FFA를 RS/RE로 구분하고 하위 영역인 VC의 RS/RE와의 연결성을 분석하였다.


t7 <- read.csv("data/bg_summary_gc_nsm.csv", header = T)
t7_l <- t7
t7_l$subj = factor(t7_l$subj, levels=c("primeRecon_2", "primeRecon_4",
                                       "primeRecon_5", "primeRecon_6",
                                       "primeRecon_7", "primeRecon_8",
                                       "primeRecon_9", "primeRecon_10",
                                       "primeRecon_11", "primeRecon_13",
                                       "primeRecon_14", "primeRecon_15",
                                       "primeRecon_16", "primeRecon_17",
                                       "primeRecon_18", "primeRecon_19",
                                       "primeRecon_20", "primeRecon_21",
                                       "primeRecon_22", "primeRecon_23",
                                       "primeRecon_24"),
                   labels=c(2,4,5,6,7,8,9,10,11,13,14,15,16,17,18,19,20,21,22,23,24))
t7_l <- t7_l %>% filter(subj!=99)
t7_l$ffa_type = factor(t7_l$ffa_type, levels = c("ffa_rs","ffa_re"))
t7_l$v_type = factor(t7_l$v_type, levels = c("v1_rs","v2_rs","v3_rs","v4_rs","all_v_rs",
                                             "v1_re","v2_re","v3_re","v4_re","all_v_re"),
                     labels=c("v1_rs","v2_rs","v3_rs","v4_rs","all_rs",
                              "v1_re","v2_re","v3_re","v4_re","all_re"))
t7_l <- t7_l %>% separate(col=v_type, sep="_", into=c("vRoi", "rsCon"))
t7_l$vRoi <- factor(t7_l$vRoi, levels = c("v1", "v2", "v3", "v4", "all"),
                    labels=c("v1", "v2", "v3", "v4", "all_vc"))
t7_l$rsCon <- factor(t7_l$rsCon, levels = c("rs", "re"))

glimpse(t7_l, width = 70)
## Rows: 400
## Columns: 14
## $ subj       <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ ffa_type   <fct> ffa_rs, ffa_rs, ffa_rs, ffa_rs, ffa_rs, ffa_rs, f…
## $ vRoi       <fct> v1, v2, v3, v4, all_vc, v1, v2, v3, v4, all_vc, v…
## $ rsCon      <fct> rs, rs, rs, rs, rs, re, re, re, re, re, rs, rs, r…
## $ Zr         <dbl> 0.4391581, 0.4276419, 0.4461436, 0.4577337, 0.499…
## $ pval       <dbl> 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.00000…
## $ vc2ffa_F   <dbl> 53.3050898, 36.2265324, 43.4217104, 66.5933437, 5…
## $ vc2ffa_Fp  <dbl> 3.597123e-13, 1.958553e-09, 5.151246e-11, 4.44089…
## $ vc2ffa_lag <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ vc2ffa_cv  <dbl> 15.17547, 15.17547, 15.17547, 15.17547, 15.17547,…
## $ ffa2vc_F   <dbl> 0.7835378, 36.2771629, 21.8846564, 13.5880824, 5.…
## $ ffa2vc_Fp  <dbl> 3.761288e-01, 1.908890e-09, 3.017064e-06, 2.31455…
## $ ffa2vc_lag <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ ffa2vc_cv  <dbl> 15.17547, 15.17547, 15.17547, 15.17547, 15.17547,…


요약치는 다음과 같다.


# subject-level, long format
t7_allL <- t7_l %>% group_by(subj, ffa_type, vRoi, rsCon) %>%
  dplyr::summarise(Zr=mean(Zr)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'ffa_type', 'vRoi'. You can override using the `.groups` argument.
# t7_allL %>% kable(digits=2)

# subject-level, wide format
t7_allW <- t7_allL %>% spread(key=ffa_type, value = Zr)
# t7_allW %>% kable(digits=2)

t7_allW1 <- t7_allL %>% spread(key=vRoi, value = Zr)
# t7_allW1 %>% kable(digits=2)

# summary table: grand mean
t7_allG <- t7_allL %>% group_by(ffa_type, vRoi, rsCon) %>%
  summarise(Zr.m = mean(Zr), Zr.sd = sd(Zr)) %>%
  ungroup()
## `summarise()` has grouped output by 'ffa_type', 'vRoi'. You can override using the `.groups` argument.
t7_allG$Zr.se <- Rmisc::summarySEwithin(data = t7_allL, measurevar = "Zr", 
                                        idvar = "subj", withinvars = c("vRoi", "rsCon", "ffa_type"))$se
t7_allG$Zr.ci <- Rmisc::summarySEwithin(data = t7_allL, measurevar = "Zr", 
                                        idvar = "subj", withinvars = c("vRoi", "rsCon", "ffa_type"))$ci
t7_allG <- t7_allG %>% 
  mutate(lower.ci = Zr.m-Zr.ci,
         upper.ci = Zr.m+Zr.ci,
         lower.se = Zr.m-Zr.se,
         upper.se = Zr.m+Zr.se)
t7_allG<-t7_allG %>% select(vRoi, ffa_type, rsCon, Zr.m, Zr.sd, Zr.se, Zr.ci, lower.ci, upper.ci, lower.se, upper.se)
t7_allG %>% select(vRoi, ffa_type, rsCon, Zr.m) %>% spread(ffa_type, Zr.m) %>% kable(digits=2)
vRoi rsCon ffa_rs ffa_re
v1 rs 0.32 0.22
v1 re 0.23 0.20
v2 rs 0.39 0.26
v2 re 0.28 0.25
v3 rs 0.44 0.28
v3 re 0.32 0.27
v4 rs 0.44 0.29
v4 re 0.36 0.28
all_vc rs 0.43 0.29
all_vc re 0.34 0.29


전반적으로 높은 연결성을 보인다. 이렇게 높은 연결성은 무엇에 의해 나타나는가? prediction signal?, granger causality analysis를 해볼 수 있을까?


어쨌든, 가설과 다르게, FFA의 RS와 VC의 RS 간 연결성이 FFA RS - VC RE 연결성보다 강하게 나타났다. 흥미롭게도, v3과 v4에서는 FFA RS - VC RE 간 연결성이 FFA RE - VC RE 보다 높게 나타났다. 이는 PC 모형에서 Error unit과 Rerepresentation Unit 간의 연결로 해석할 수 있는 가능성을 보여준다. RS가 Representation Unit (RU)이라고 한다면, FFA의 RU와 VC의 RU 간 강한 연결성은 자극에 대한 Representation 자체, 또는 sensory forward signal에 기인할 수 있다 (PC 모델 그림에서는 이 연결을 표시하지 않고 있는데, Forward Sensory Signal이 있는 것으로 보는게 타당할 것이다.). 아니면.. Prediction의 Feedback Signal이거나… RE를 Error Unit (EU)이라고 한다면, FFA의 RU와 VC의 EU 간 강한 연결성은 Prediction Error에 의한 Forward Signal일 수 있다. 이 가설을 검증하기 위해서는 VC RU <-> FFA RU, VC EU -> FFA EU 의 Effective Connectivity 패턴이 나타날 가능성이 있다 (VC RU <-> FFA RU 에서의 방향성은 추가적인 검토가 필요하다. Forward Signal인가? Feedback Signal인가?, 또 V1 & V2에서는 왜 안나타나는가? 너무 먼 hierarchy라서?). 일단 어느 방향으로 가정을 하든, Granger Causality Analysis를 통해 방향성이 존재하는지 살펴볼 수는 있을 것이다.


또 한가지 흥미로운 점은 FFA의 RS 복셀들은 전반적으로 하위 영역의 RS/RE와 높은 연결성을 보이는데, FFA의 RE 복셀들은 조건간 차이를 보이지 않는다는 점이다. Error unit이 forward signal을 투사하는 방식으로 상위 영역과 연결되어 있다는 점을 적용하여 설명할 수 있는 것으로 보인다.


한가지 고려해야하는 점은, Results 1의 분석에서 RE의 효과 크기나 RE 복셀의 비율 자체는 상위 영역과의 연결성이 강할수록 감소했다는 점이다. 물론 여기서 상위 영역인 FFA는 RS/RE를 구분하지 않았다. RS 복셀이 RE 복셀보다 더 많다는 것을 생각해보면 RS 복셀들의 time course가 전체 FFA time course에 반영될 가능성이 크다. 추가 분석으로 FFA RS, FFA RE를 기준으로 Results 1의 분석을 다시 해볼 필요가 있다.


t7.all.plot7.1 <- ggplot(data=t7_allL, aes(x=ffa_type, y=Zr, fill=rsCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_point(data=t7_allL, aes(x=ffa_type, y=Zr, fill=rsCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t7_allG, aes(x = ffa_type, y=Zr.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~vRoi, scales="free_x", space = "free",
             labeller = labeller(vRoi = c("v1" = "V1","v2" = "V2","v3" = "V3","v4" = "V4", "all_vc" = "all VC"))) +
  scale_x_discrete(labels=c("RS","RE")) +
  scale_fill_manual(values = c("#feb24c", "#91bfdb"), # c("#feb24c", "#91bfdb"),
                    labels = c("VC RS voxel", "VC RE voxel")) +
  coord_cartesian(ylim = c(-0.2, 0.7), clip = "on") +
  labs(x = "FFA", y = "Background Connectivity (Zr)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# t7.all.plot7.1
#ggsave("bg_ffa_plot1.jpg", plot = t7.all.plot7.1, width=12, height=8, unit='in', dpi=600)

t7.all.plot7.2 <- ggplot(data=t7_allL, aes(x=vRoi, y=Zr, fill=ffa_type)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_point(data=t7_allL, aes(x=vRoi, y=Zr, fill=ffa_type), position = position_dodge(width=0.8),
           size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t7_allG, aes(x =vRoi, y=Zr.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~rsCon, scales="free_x", space = "free",
             labeller = labeller(rsCon = c("rs" = "Visual RS voxel","re" = "Visual RE voxel"))) +
  scale_x_discrete(labels=c("V1","V3","V3","V4","all VC")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),  # , "#235796"
                    labels = c("FFA RS voxel", "FFA RE voxel")) +
  coord_cartesian(ylim = c(-0.2, 0.7), clip = "on") +
  labs(x = "Visucal Cortex", y = "Background Connectivity (Zr)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# t7.all.plot7.2
#ggsave("bg_ffa_plot2.jpg", plot = t7.all.plot7.2, width=12, height=8, unit='in', dpi=600)

t7.all.plot7.3 <- ggplot(data=t7_allL, aes(x=rsCon, y=Zr, fill=ffa_type)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_point(data=t7_allL, aes(x=rsCon, y=Zr, fill=ffa_type), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t7_allG, aes(x =rsCon, y=Zr.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~vRoi, scales="free_x", space = "free",
             labeller = labeller(vRoi = c("v1" = "V1","v2" = "V2","v3" = "V3","v4" = "V4", "all_vc" = "all VC"))) +
  scale_x_discrete(labels=c("RS","RE")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),  # , "#235796"
                    labels = c("FFA RS voxel", "FFA RE voxel")) +
  coord_cartesian(ylim = c(-0.2, 0.7), clip = "on") +
  labs(x = "Visual Cortex", y = "Background Connectivity (Zr)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
# t7.all.plot7.3
#ggsave("bg_ffa_plot3.jpg", plot = t7.all.plot7.3, width=12, height=8, unit='in', dpi=600)

# 
ggarrange(t7.all.plot7.1, t7.all.plot7.3,
          ncol=2, labels=c("A","B"))


시각 관심 영역(V1 ~ V4, all VC), FFA의 RS/RE, VC의 RS/RE를 요인으로 5 x 2 x 2 혼합 요인 분산분석을 수행하였다. 다른 것도 있지만, 주목할만한 점으로는 VC의 RS/RE와 FFA의 RS/RE 간의 상호작용이 유의하였다.


t7_allL.1 <- t7_allL
t7.aov1 <- aov_ez(id="subj", dv="Zr", data = t7_allL.1, within = c("vRoi", "rsCon", "ffa_type"))
nice(t7.aov1, es="pes") %>% kable(digits=2)
Effect df MSE F pes p.value
vRoi 2.02, 38.35 0.02 15.82 *** .454 <.001
rsCon 1, 19 0.04 7.38 * .280 .014
ffa_type 1, 19 0.06 13.27 ** .411 .002
vRoi:rsCon 2.59, 49.20 0.00 0.57 .029 .614
vRoi:ffa_type 1.45, 27.48 0.01 3.24 + .146 .069
rsCon:ffa_type 1, 19 0.01 23.73 *** .555 <.001
vRoi:rsCon:ffa_type 2.17, 41.27 0.00 2.00 .095 .144


사후 분석 결과, FFA RS 에서 VC RS 와 VC RE 간 차이가 유의하였다.구체적으로 FFA RS-VC RS의 연결성 강도가 크게 나타났다.


t7.aov1.m2 <- emmeans(t7.aov1, pairwise ~ rsCon | ffa_type + vRoi , adjust = "tukey") # adjust="bon" , "tukey", "none"
t7.aov1.m2$contrasts %>% kable(digits=3)
contrast ffa_type vRoi estimate SE df t.ratio p.value
rs - re ffa_rs v1 0.091 0.028 19 3.231 0.004
rs - re ffa_re v1 0.016 0.022 19 0.722 0.479
rs - re ffa_rs v2 0.107 0.032 19 3.282 0.004
rs - re ffa_re v2 0.007 0.022 19 0.312 0.759
rs - re ffa_rs v3 0.116 0.024 19 4.950 0.000
rs - re ffa_re v3 0.014 0.022 19 0.640 0.530
rs - re ffa_rs v4 0.080 0.021 19 3.801 0.001
rs - re ffa_re v4 0.007 0.016 19 0.416 0.682
rs - re ffa_rs all_vc 0.095 0.029 19 3.270 0.004
rs - re ffa_re all_vc -0.004 0.022 19 -0.194 0.848
# plot(t7.aov1.m2, horizontal = T, comparison=T)


V4의 VC RE 에서 FFA RS와 FFA RE 간 차이가 유의하였다. V3에서는 Marginal한 차이를 보였다. 구체적으로 VC RE-FFA RS의 연결성 강도가 VC RE-FFA RE 보다 크게 나타났다. V1, V2에서는 이런 효과가 나타나지 않았는데, 효과가 충분하지 않은 것인지, Error signal의 연결이 V1 -> V2 -> V3 -> V4로 흐르기 때문에 hierarchy 상 거리가 멀어서 그렇게 나타난 것인지 알 수 없다. 또한 이 데이터에는 Meridian Mapping이 없어서 VC ROI가 부정확할 가능성도 있다.


t7.aov1.m1 <- emmeans(t7.aov1, pairwise ~ ffa_type | rsCon + vRoi, adjust = "tukey") # adjust="bon" , "tukey", "none"
t7.aov1.m1$contrasts %>% kable(digits=3)
contrast rsCon vRoi estimate SE df t.ratio p.value
ffa_rs - ffa_re rs v1 0.099 0.031 19 3.145 0.005
ffa_rs - ffa_re re v1 0.024 0.032 19 0.763 0.455
ffa_rs - ffa_re rs v2 0.130 0.027 19 4.831 0.000
ffa_rs - ffa_re re v2 0.031 0.032 19 0.963 0.348
ffa_rs - ffa_re rs v3 0.152 0.024 19 6.233 0.000
ffa_rs - ffa_re re v3 0.050 0.026 19 1.925 0.069
ffa_rs - ffa_re rs v4 0.152 0.026 19 5.835 0.000
ffa_rs - ffa_re re v4 0.078 0.024 19 3.294 0.004
ffa_rs - ffa_re rs all_vc 0.145 0.029 19 4.962 0.000
ffa_rs - ffa_re re all_vc 0.046 0.034 19 1.329 0.200
# plot(t7.aov1.m1, horizontal = T, comparison=T)

3.2 Visual RS/RE와 상위 영역 RS/RE 간의 유효 연결성



# subject-level, long format
t8_l <- t7_l
thr = t8_l$vc2ffa_cv[1]

t8_l <- gather(t8_l, dirt, F, c(vc2ffa_F, ffa2vc_F), factor_key=TRUE)
t8_l$dirt <- factor(t8_l$dirt, levels = c('vc2ffa_F', 'ffa2vc_F'), labels=c('forward','back'))

t8_allL <- t8_l %>% group_by(subj, ffa_type, vRoi, rsCon, dirt) %>%
  dplyr::summarise(F=mean(F)) %>%
  ungroup()
## `summarise()` has grouped output by 'subj', 'ffa_type', 'vRoi', 'rsCon'. You can override using the `.groups` argument.
# t8_allL %>% kable(digits=2)

# subject-level, wide format
t8_allW <- t8_allL %>% spread(key=ffa_type, value = F)
# t8_allW %>% kable(digits=2)

t8_allW1 <- t8_allL %>% spread(key=vRoi, value = F)
# t8_allW1 %>% kable(digits=2)

# t8_allW2 <- t8_allL %>% spread(key=dirt, value = F)
# t8_allW2 %>% kable(digits=2)

# summary table: grand mean
t8_allG <- t8_allL %>% group_by(ffa_type, vRoi, rsCon, dirt) %>%
  summarise(F.m = mean(F), F.sd = sd(F)) %>%
  ungroup()
## `summarise()` has grouped output by 'ffa_type', 'vRoi', 'rsCon'. You can override using the `.groups` argument.
t8_allG$F.se <- Rmisc::summarySEwithin(data = t8_allL, measurevar = "F", 
                                        idvar = "subj", withinvars = c("vRoi", "rsCon", "ffa_type", "dirt"))$se
t8_allG$F.ci <- Rmisc::summarySEwithin(data = t8_allL, measurevar = "F", 
                                        idvar = "subj", withinvars = c("vRoi", "rsCon", "ffa_type", "dirt"))$ci
t8_allG <- t8_allG %>% 
  mutate(lower.ci = F.m-F.ci,
         upper.ci = F.m+F.ci,
         lower.se = F.m-F.se,
         upper.se = F.m+F.se)
t8_allG<-t8_allG %>% select(vRoi, ffa_type, rsCon, dirt, F.m, F.sd, F.se, F.ci, lower.ci, upper.ci, lower.se, upper.se)
t8_allG %>% select(vRoi, ffa_type, rsCon, dirt, F.m) %>% spread(ffa_type, F.m) %>% kable(digits=2)
vRoi rsCon dirt ffa_rs ffa_re
v1 rs forward 32.08 27.32
v1 rs back 35.31 22.35
v1 re forward 19.30 18.54
v1 re back 28.95 20.94
v2 rs forward 47.26 39.40
v2 rs back 40.94 24.37
v2 re forward 34.73 29.55
v2 re back 30.03 24.48
v3 rs forward 56.44 50.85
v3 rs back 45.42 26.85
v3 re forward 51.78 50.30
v3 re back 35.61 24.25
v4 rs forward 64.02 52.95
v4 rs back 41.97 22.26
v4 re forward 52.55 43.57
v4 re back 52.37 28.86
all_vc rs forward 55.04 46.55
all_vc rs back 37.50 22.96
all_vc re forward 47.40 41.95
all_vc re back 34.44 23.61


t8_allL.1 <- t8_allL %>% filter(ffa_type == "ffa_rs")
t8_allG.1 <- t8_allG %>% filter(ffa_type == "ffa_rs")
t8.all.plot8.1 <- ggplot(data=t8_allL.1, aes(x=dirt, y=F, fill=rsCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
  geom_point(data=t8_allL.1, aes(x=dirt, y=F, fill=rsCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t8_allG.1, aes(x = dirt, y=F.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~vRoi, scales="free_x", space = "free",
             labeller = labeller(vRoi = c("v1" = "V1","v2" = "V2","v3" = "V3","v4" = "V4", "all_vc" = "all VC"))) +
  scale_x_discrete(labels=c("to FFA","from FFA")) +
  scale_fill_manual(values = c("#feb24c", "#91bfdb"), # c("#feb24c", "#91bfdb"),
                    labels = c("VC RS voxel", "VC RE voxel")) +
  coord_cartesian(ylim = c(0, 100), clip = "on") +
  labs(x = "FFA RS", y = "Granger Causality F Value") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
#t8.all.plot8.1

t8_allL.2 <- t8_allL %>% filter(ffa_type == "ffa_re")
t8_allG.2 <- t8_allG %>% filter(ffa_type == "ffa_re")
t8.all.plot8.2 <- ggplot(data=t8_allL.2, aes(x=dirt, y=F, fill=rsCon)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
  geom_point(data=t8_allL.2, aes(x=dirt, y=F, fill=rsCon), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t8_allG.2, aes(x = dirt, y=F.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~vRoi, scales="free_x", space = "free",
             labeller = labeller(vRoi = c("v1" = "V1","v2" = "V2","v3" = "V3","v4" = "V4", "all_vc" = "all VC"))) +
  scale_x_discrete(labels=c("to FFA","from FFA")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"), # c("#feb24c", "#91bfdb"),
                    labels = c("VC RS voxel", "VC RE voxel")) +
  coord_cartesian(ylim = c(0, 100), clip = "on") +
  labs(x = "FFA RE", y = "Granger Causality F Value") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
#t8.all.plot8.2
ggarrange(t8.all.plot8.1, t8.all.plot8.2,
          ncol=2, labels=c("A","B"))


t8_allL.1 <- t8_allL %>% filter(ffa_type == "ffa_rs")
t8_allG.1 <- t8_allG %>% filter(ffa_type == "ffa_rs")
t8.all.plot8.3 <- ggplot(data=t8_allL.1, aes(x=rsCon, y=F, fill=dirt)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
  geom_point(data=t8_allL.1, aes(x=rsCon, y=F, fill=dirt), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t8_allG.1, aes(x = rsCon, y=F.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~vRoi, scales="free_x", space = "free",
             labeller = labeller(vRoi = c("v1" = "V1","v2" = "V2","v3" = "V3","v4" = "V4", "all_vc" = "all VC"))) +
  scale_x_discrete(labels=c("RS","RE")) +
  scale_fill_manual(values = c("#feb24c", "#91bfdb"), # c("#feb24c", "#91bfdb"),
                    labels = c("Forward to FFA RS", "Back from FFA RS")) +
  coord_cartesian(ylim = c(0, 100), clip = "on") +
  labs(x = "VC RS/RE", y = "Granger Causality F Value") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
#t8.all.plot8.3

t8_allL.2 <- t8_allL %>% filter(ffa_type == "ffa_re")
t8_allG.2 <- t8_allG %>% filter(ffa_type == "ffa_re")
t8.all.plot8.4 <- ggplot(data=t8_allL.2, aes(x=rsCon, y=F, fill=dirt)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
  geom_point(data=t8_allL.2, aes(x=rsCon, y=F, fill=dirt), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t8_allG.2, aes(x = rsCon, y=F.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~vRoi, scales="free_x", space = "free",
             labeller = labeller(vRoi = c("v1" = "V1","v2" = "V2","v3" = "V3","v4" = "V4", "all_vc" = "all VC"))) +
  scale_x_discrete(labels=c("RS","RE")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"), # c("#feb24c", "#91bfdb"),
                    labels = c("Forward to FFA RE", "Back from FFA RE")) +
  coord_cartesian(ylim = c(0, 100), clip = "on") +
  labs(x = "VC RS/RE", y = "Granger Causality F Value") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
#t8.all.plot8.4
ggarrange(t8.all.plot8.3, t8.all.plot8.4,
          ncol=2, labels=c("A","B"))


# t8_allL.3 <- t8_allL %>% filter(vRoi == "v1")
# t8_allG.3 <- t8_allG %>% filter(vRoi == "v1")
# t8.all.plot8.5 <- ggplot(data=t8_allL.3, aes(x=rsCon, y=F, fill=dirt)) +
#   stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
#   geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
#   geom_point(data=t8_allL.3, aes(x=rsCon, y=F, fill=dirt), position = position_dodge(width=0.8),
#              size=2, show.legend=F, color="gray80") +
#   geom_pointrange(data=t8_allG.3, aes(x = rsCon, y=F.m, ymin = lower.se, ymax = upper.se),
#                   position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
#   facet_grid(.~ffa_type, scales="free_x", space = "free",
#              labeller = labeller(ffa_type = c("ffa_rs" = "FFA RS","ffa_re" = "FFA RE"))) +
#   scale_x_discrete(labels=c("RS","RE")) +
#   scale_fill_manual(values = c("#feb24c", "#91bfdb"), # c("#feb24c", "#91bfdb"),
#                     labels = c("FeedForward", "FeedBack")) +
#   coord_cartesian(ylim = c(0, 100), clip = "on") +
#   labs(x = "V3", y = "Granger Causality F Value") +
#   theme_bw(base_size = 18) +
#   theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
#         axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
#         axis.line=element_line(),
#         strip.text.x = element_text(face = "plain", size = 15, color = "black"),
#         strip.background = element_blank(),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.border = element_blank(),
#         panel.spacing=unit(1, "lines"),
#         plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
#         legend.title = element_blank(),
#         legend.position="top")
# #t8.all.plot8.5
# 
# t8.all.plot8.6 <- ggplot(data=t8_allL.3, aes(x=ffa_type, y=F, fill=dirt)) +
#   stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
#   geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
#   geom_point(data=t8_allL.3, aes(x=ffa_type, y=F, fill=dirt), position = position_dodge(width=0.8),
#              size=2, show.legend=F, color="gray80") +
#   geom_pointrange(data=t8_allG.3, aes(x = ffa_type, y=F.m, ymin = lower.se, ymax = upper.se),
#                   position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
#   facet_grid(.~rsCon, scales="free_x", space = "free",
#              labeller = labeller(rsCon = c("rs" = "V1 RS","re" = "V1 RE"))) +
#   scale_x_discrete(labels=c("RS","RE")) +
#   scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"), # c("#feb24c", "#91bfdb"),
#                     labels = c("FeedForward", "Feedback")) +
#   coord_cartesian(ylim = c(0, 120), clip = "on") +
#   labs(x = "FFA", y = "Granger Causality F Value") +
#   theme_bw(base_size = 18) +
#   theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
#         axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
#         axis.line=element_line(),
#         strip.text.x = element_text(face = "plain", size = 15, color = "black"),
#         strip.background = element_blank(),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.border = element_blank(),
#         panel.spacing=unit(1, "lines"),
#         plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
#         legend.title = element_blank(),
#         legend.position="top")
# #t8.all.plot8.6
# ggarrange(t8.all.plot8.5, t8.all.plot8.6,
#           ncol=2, labels=c("A","B"))


# t8_allL.3 <- t8_allL %>% filter(vRoi == "v2")
# t8_allG.3 <- t8_allG %>% filter(vRoi == "v2")
# t8.all.plot8.5 <- ggplot(data=t8_allL.3, aes(x=rsCon, y=F, fill=dirt)) +
#   stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
#   geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
#   geom_point(data=t8_allL.3, aes(x=rsCon, y=F, fill=dirt), position = position_dodge(width=0.8),
#              size=2, show.legend=F, color="gray80") +
#   geom_pointrange(data=t8_allG.3, aes(x = rsCon, y=F.m, ymin = lower.se, ymax = upper.se),
#                   position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
#   facet_grid(.~ffa_type, scales="free_x", space = "free",
#              labeller = labeller(ffa_type = c("ffa_rs" = "FFA RS","ffa_re" = "FFA RE"))) +
#   scale_x_discrete(labels=c("RS","RE")) +
#   scale_fill_manual(values = c("#feb24c", "#91bfdb"), # c("#feb24c", "#91bfdb"),
#                     labels = c("FeedForward", "FeedBack")) +
#   coord_cartesian(ylim = c(0, 100), clip = "on") +
#   labs(x = "V2", y = "Granger Causality F Value") +
#   theme_bw(base_size = 18) +
#   theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
#         axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
#         axis.line=element_line(),
#         strip.text.x = element_text(face = "plain", size = 15, color = "black"),
#         strip.background = element_blank(),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.border = element_blank(),
#         panel.spacing=unit(1, "lines"),
#         plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
#         legend.title = element_blank(),
#         legend.position="top")
# #t8.all.plot8.5
# 
# t8.all.plot8.6 <- ggplot(data=t8_allL.3, aes(x=ffa_type, y=F, fill=dirt)) +
#   stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
#   geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
#   geom_point(data=t8_allL.3, aes(x=ffa_type, y=F, fill=dirt), position = position_dodge(width=0.8),
#              size=2, show.legend=F, color="gray80") +
#   geom_pointrange(data=t8_allG.3, aes(x = ffa_type, y=F.m, ymin = lower.se, ymax = upper.se),
#                   position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
#   facet_grid(.~rsCon, scales="free_x", space = "free",
#              labeller = labeller(rsCon = c("rs" = "V2 RS","re" = "V2 RE"))) +
#   scale_x_discrete(labels=c("RS","RE")) +
#   scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"), # c("#feb24c", "#91bfdb"),
#                     labels = c("FeedForward", "Feedback")) +
#   coord_cartesian(ylim = c(0, 120), clip = "on") +
#   labs(x = "FFA", y = "Granger Causality F Value") +
#   theme_bw(base_size = 18) +
#   theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
#         axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
#         axis.line=element_line(),
#         strip.text.x = element_text(face = "plain", size = 15, color = "black"),
#         strip.background = element_blank(),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.border = element_blank(),
#         panel.spacing=unit(1, "lines"),
#         plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
#         legend.title = element_blank(),
#         legend.position="top")
# #t8.all.plot8.6
# ggarrange(t8.all.plot8.5, t8.all.plot8.6,
#           ncol=2, labels=c("A","B"))


# t8_allL.3 <- t8_allL %>% filter(vRoi == "v3")
# t8_allG.3 <- t8_allG %>% filter(vRoi == "v3")
# t8.all.plot8.5 <- ggplot(data=t8_allL.3, aes(x=rsCon, y=F, fill=dirt)) +
#   stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
#   geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
#   geom_point(data=t8_allL.3, aes(x=rsCon, y=F, fill=dirt), position = position_dodge(width=0.8),
#              size=2, show.legend=F, color="gray80") +
#   geom_pointrange(data=t8_allG.3, aes(x = rsCon, y=F.m, ymin = lower.se, ymax = upper.se),
#                   position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
#   facet_grid(.~ffa_type, scales="free_x", space = "free",
#              labeller = labeller(ffa_type = c("ffa_rs" = "FFA RS","ffa_re" = "FFA RE"))) +
#   scale_x_discrete(labels=c("RS","RE")) +
#   scale_fill_manual(values = c("#feb24c", "#91bfdb"), # c("#feb24c", "#91bfdb"),
#                     labels = c("FeedForward", "FeedBack")) +
#   coord_cartesian(ylim = c(0, 100), clip = "on") +
#   labs(x = "V3", y = "Granger Causality F Value") +
#   theme_bw(base_size = 18) +
#   theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
#         axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
#         axis.line=element_line(),
#         strip.text.x = element_text(face = "plain", size = 15, color = "black"),
#         strip.background = element_blank(),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.border = element_blank(),
#         panel.spacing=unit(1, "lines"),
#         plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
#         legend.title = element_blank(),
#         legend.position="top")
# #t8.all.plot8.5
# 
# t8.all.plot8.6 <- ggplot(data=t8_allL.3, aes(x=ffa_type, y=F, fill=dirt)) +
#   stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
#   geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
#   geom_point(data=t8_allL.3, aes(x=ffa_type, y=F, fill=dirt), position = position_dodge(width=0.8),
#              size=2, show.legend=F, color="gray80") +
#   geom_pointrange(data=t8_allG.3, aes(x = ffa_type, y=F.m, ymin = lower.se, ymax = upper.se),
#                   position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
#   facet_grid(.~rsCon, scales="free_x", space = "free",
#              labeller = labeller(rsCon = c("rs" = "V3 RS","re" = "V3 RE"))) +
#   scale_x_discrete(labels=c("RS","RE")) +
#   scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"), # c("#feb24c", "#91bfdb"),
#                     labels = c("FeedForward", "Feedback")) +
#   coord_cartesian(ylim = c(0, 120), clip = "on") +
#   labs(x = "FFA", y = "Granger Causality F Value") +
#   theme_bw(base_size = 18) +
#   theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
#         axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
#         axis.line=element_line(),
#         strip.text.x = element_text(face = "plain", size = 15, color = "black"),
#         strip.background = element_blank(),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.border = element_blank(),
#         panel.spacing=unit(1, "lines"),
#         plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
#         legend.title = element_blank(),
#         legend.position="top")
# #t8.all.plot8.6
# ggarrange(t8.all.plot8.5, t8.all.plot8.6,
#           ncol=2, labels=c("A","B"))


t8_allL.3 <- t8_allL %>% filter(vRoi == "v4")
t8_allG.3 <- t8_allG %>% filter(vRoi == "v4")
t8.all.plot8.5 <- ggplot(data=t8_allL.3, aes(x=rsCon, y=F, fill=dirt)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
  geom_point(data=t8_allL.3, aes(x=rsCon, y=F, fill=dirt), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t8_allG.3, aes(x = rsCon, y=F.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~ffa_type, scales="free_x", space = "free",
             labeller = labeller(ffa_type = c("ffa_rs" = "FFA RS","ffa_re" = "FFA RE"))) +
  scale_x_discrete(labels=c("RS","RE")) +
  scale_fill_manual(values = c("#feb24c", "#91bfdb"), # c("#feb24c", "#91bfdb"),
                    labels = c("FeedForward", "FeedBack")) +
  coord_cartesian(ylim = c(0, 100), clip = "on") +
  labs(x = "V4", y = "Granger Causality F Value") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
#t8.all.plot8.5

t8.all.plot8.6 <- ggplot(data=t8_allL.3, aes(x=ffa_type, y=F, fill=dirt)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
  geom_point(data=t8_allL.3, aes(x=ffa_type, y=F, fill=dirt), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t8_allG.3, aes(x = ffa_type, y=F.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~rsCon, scales="free_x", space = "free",
             labeller = labeller(rsCon = c("rs" = "V4 RS","re" = "V4 RE"))) +
  scale_x_discrete(labels=c("RS","RE")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"), # c("#feb24c", "#91bfdb"),
                    labels = c("FeedForward", "Feedback")) +
  coord_cartesian(ylim = c(0, 120), clip = "on") +
  labs(x = "FFA", y = "Granger Causality F Value") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
#t8.all.plot8.6
ggarrange(t8.all.plot8.5, t8.all.plot8.6,
          ncol=2, labels=c("A","B"))

t8.all.plot8.7 <- ggplot(data=t8_allL.3, aes(x=rsCon, y=F, fill=ffa_type)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
  geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
  geom_point(data=t8_allL.3, aes(x=rsCon, y=F, fill=ffa_type), position = position_dodge(width=0.8),
             size=2, show.legend=F, color="gray80") +
  geom_pointrange(data=t8_allG.3, aes(x = rsCon, y=F.m, ymin = lower.se, ymax = upper.se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  facet_grid(.~dirt, scales="free_x", space = "free",
             labeller = labeller(dirt = c("forward" = "V4 to FFA","back" = "FFA to V4"))) +
  scale_x_discrete(labels=c("RS","RE")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"), # c("#feb24c", "#91bfdb"),
                    labels = c("FFA RS", "FFA RE")) +
  coord_cartesian(ylim = c(0, 120), clip = "on") +
  labs(x = "V4", y = "Granger Causality F Value") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
        legend.title = element_blank(),
        legend.position="top")
t8.all.plot8.7


# t8_allL.3 <- t8_allL %>% filter(vRoi == "all_vc")
# t8_allG.3 <- t8_allG %>% filter(vRoi == "all_vc")
# t8.all.plot8.5 <- ggplot(data=t8_allL.3, aes(x=rsCon, y=F, fill=dirt)) +
#   stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
#   geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
#   geom_point(data=t8_allL.3, aes(x=rsCon, y=F, fill=dirt), position = position_dodge(width=0.8),
#              size=2, show.legend=F, color="gray80") +
#   geom_pointrange(data=t8_allG.3, aes(x = rsCon, y=F.m, ymin = lower.se, ymax = upper.se),
#                   position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
#   facet_grid(.~ffa_type, scales="free_x", space = "free",
#              labeller = labeller(ffa_type = c("ffa_rs" = "FFA RS","ffa_re" = "FFA RE"))) +
#   scale_x_discrete(labels=c("RS","RE")) +
#   scale_fill_manual(values = c("#feb24c", "#91bfdb"), # c("#feb24c", "#91bfdb"),
#                     labels = c("FeedForward", "FeedBack")) +
#   coord_cartesian(ylim = c(0, 100), clip = "on") +
#   labs(x = "all VC", y = "Granger Causality F Value") +
#   theme_bw(base_size = 18) +
#   theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
#         axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
#         axis.line=element_line(),
#         strip.text.x = element_text(face = "plain", size = 15, color = "black"),
#         strip.background = element_blank(),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.border = element_blank(),
#         panel.spacing=unit(1, "lines"),
#         plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
#         legend.title = element_blank(),
#         legend.position="top")
# #t8.all.plot8.5
# 
# t8.all.plot8.6 <- ggplot(data=t8_allL.3, aes(x=ffa_type, y=F, fill=dirt)) +
#   stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, width = 0.8,  size = 1.02) +
#   geom_hline(yintercept=thr, linetype='dashed', color='gray60', alpha =1, size=1) +
#   geom_point(data=t8_allL.3, aes(x=ffa_type, y=F, fill=dirt), position = position_dodge(width=0.8),
#              size=2, show.legend=F, color="gray80") +
#   geom_pointrange(data=t8_allG.3, aes(x = ffa_type, y=F.m, ymin = lower.se, ymax = upper.se),
#                   position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
#   facet_grid(.~rsCon, scales="free_x", space = "free",
#              labeller = labeller(rsCon = c("rs" = "all VC RS","re" = "all VC RE"))) +
#   scale_x_discrete(labels=c("RS","RE")) +
#   scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"), # c("#feb24c", "#91bfdb"),
#                     labels = c("FeedForward", "Feedback")) +
#   coord_cartesian(ylim = c(0, 120), clip = "on") +
#   labs(x = "FFA", y = "Granger Causality F Value") +
#   theme_bw(base_size = 18) +
#   theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
#         axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
#         axis.line=element_line(),
#         strip.text.x = element_text(face = "plain", size = 15, color = "black"),
#         strip.background = element_blank(),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.border = element_blank(),
#         panel.spacing=unit(1, "lines"),
#         plot.margin = margin(1, 0.3, 1, 0.3, "cm"), 
#         legend.title = element_blank(),
#         legend.position="top")
# #t8.all.plot8.6
# ggarrange(t8.all.plot8.5, t8.all.plot8.6,
#           ncol=2, labels=c("A","B"))


시각 관심 영역(V4, all VC), FFA의 RS/RE, VC의 RS/RE를 요인으로 5 x 2 x 2 혼합 요인 분산분석을 수행하였다. 다른 것도 있지만, 주목할만한 점으로는 VC의 RS/RE와 FFA의 RS/RE 간의 상호작용이 유의하였다.


t8_allL.3 <- t8_allL %>% filter(vRoi == "v4")
t8.aov1 <- aov_ez(id="subj", dv="F", data = t8_allL.3, within = c("rsCon", "ffa_type", "dirt"))
nice(t8.aov1, es="pes") %>% kable(digits=2)
Effect df MSE F pes p.value
rsCon 1, 19 558.09 0.07 .003 .799
ffa_type 1, 19 1891.29 5.29 * .218 .033
dirt 1, 19 7861.35 1.45 .071 .243
rsCon:ffa_type 1, 19 174.29 0.04 .002 .841
rsCon:dirt 1, 19 1388.17 2.58 .120 .125
ffa_type:dirt 1, 19 1613.40 0.83 .042 .373
rsCon:ffa_type:dirt 1, 19 183.02 0.47 .024 .500


사후 분석 결과, FFA RS 에서 VC RS 와 VC RE 간 차이가 유의하였다.구체적으로 FFA RS-VC RS의 연결성 강도가 크게 나타났다.


t8.aov1.m2 <- emmeans(t8.aov1, pairwise ~ dirt | rsCon + ffa_type , adjust = "tukey") # adjust="bon" , "tukey", "none"
t8.aov1.m2$contrasts %>% kable(digits=3)
contrast rsCon ffa_type estimate SE df t.ratio p.value
forward - back rs ffa_rs 22.051 18.798 19 1.173 0.255
forward - back re ffa_rs 0.184 18.694 19 0.010 0.992
forward - back rs ffa_re 30.691 16.040 19 1.913 0.071
forward - back re ffa_re 14.709 12.019 19 1.224 0.236
# plot(t8.aov1.m2, horizontal = T, comparison=T)

t8.aov1.m2 <- emmeans(t8.aov1, pairwise ~ dirt | ffa_type  + rsCon  , adjust = "tukey") # adjust="bon" , "tukey", "none"
t8.aov1.m2$contrasts %>% kable(digits=3)
contrast ffa_type rsCon estimate SE df t.ratio p.value
forward - back ffa_rs rs 22.051 18.798 19 1.173 0.255
forward - back ffa_re rs 30.691 16.040 19 1.913 0.071
forward - back ffa_rs re 0.184 18.694 19 0.010 0.992
forward - back ffa_re re 14.709 12.019 19 1.224 0.236
# plot(t8.aov1.m2, horizontal = T, comparison=T)


t8.aov1.m2 <- emmeans(t8.aov1, pairwise ~ rsCon | dirt + ffa_type , adjust = "tukey") # adjust="bon" , "tukey", "none"
t8.aov1.m2$contrasts %>% kable(digits=3)
contrast dirt ffa_type estimate SE df t.ratio p.value
rs - re forward ffa_rs 11.472 10.617 19 1.081 0.293
rs - re back ffa_rs -10.396 5.482 19 -1.896 0.073
rs - re forward ffa_re 9.379 7.463 19 1.257 0.224
rs - re back ffa_re -6.603 5.647 19 -1.169 0.257
# plot(t8.aov1.m2, horizontal = T, comparison=T)

t8.aov1.m2 <- emmeans(t8.aov1, pairwise ~ rsCon | ffa_type + dirt , adjust = "tukey") # adjust="bon" , "tukey", "none"
t8.aov1.m2$contrasts %>% kable(digits=3)
contrast ffa_type dirt estimate SE df t.ratio p.value
rs - re ffa_rs forward 11.472 10.617 19 1.081 0.293
rs - re ffa_re forward 9.379 7.463 19 1.257 0.224
rs - re ffa_rs back -10.396 5.482 19 -1.896 0.073
rs - re ffa_re back -6.603 5.647 19 -1.169 0.257
# plot(t8.aov1.m2, horizontal = T, comparison=T)

t8.aov1.m2 <- emmeans(t8.aov1, pairwise ~ ffa_type | rsCon + dirt , adjust = "tukey") # adjust="bon" , "tukey", "none"
t8.aov1.m2$contrasts %>% kable(digits=3)
contrast rsCon dirt estimate SE df t.ratio p.value
ffa_rs - ffa_re rs forward 11.072 11.491 19 0.964 0.347
ffa_rs - ffa_re re forward 8.979 7.788 19 1.153 0.263
ffa_rs - ffa_re rs back 19.712 9.533 19 2.068 0.053
ffa_rs - ffa_re re back 23.505 10.130 19 2.320 0.032
# plot(t8.aov1.m2, horizontal = T, comparison=T)

t8.aov1.m2 <- emmeans(t8.aov1, pairwise ~ ffa_type | dirt + rsCon , adjust = "tukey") # adjust="bon" , "tukey", "none"
t8.aov1.m2$contrasts %>% kable(digits=3)
contrast dirt rsCon estimate SE df t.ratio p.value
ffa_rs - ffa_re forward rs 11.072 11.491 19 0.964 0.347
ffa_rs - ffa_re back rs 19.712 9.533 19 2.068 0.053
ffa_rs - ffa_re forward re 8.979 7.788 19 1.153 0.263
ffa_rs - ffa_re back re 23.505 10.130 19 2.320 0.032
# plot(t8.aov1.m2, horizontal = T, comparison=T)


~~